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These notes are an introduction to the theory of stochastic pro- 
cesses based on several sources. The presentation mainly follows 
the books of van Kampen [5 J and Wio [6], except for the introduc- 
tion, which is taken from the book of Gardiner [2] and the parts 
devoted to the Langevin equation and the methods for solving 
Langevin and Fokker-Planck equations, which are based on the 
book of Risken [4] . 
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1 Historical introduction 



Theoretical science up to the end of the nineteenth century can be roughly 
viewed as the study of solutions of differential equations and the modelling of 
natural phenomena by deterministic solutions of these differential equations. 
It was at that time commonly thought that if all initial (and contour) data 
could only be collected, one would be able to predict the future with certainty. 

We now know that this is not so, in at least two ways. First, the ad- 
vent of quantum mechanics gave rise to a new physics, which had as an 
essential ingredient a purely statistical element (the measurement process). 
Secondly, the concept of chaos has arisen, in which even quite simple differ- 
ential equations have the rather alarming property of giving rise to essentially 
unpredictable behaviours. 

Chaos and quantum mechanics are not the subject of these notes, but 
we shall deal with systems were limited predictability arises in the form of 
fluctuations due to the finite number of their discrete constituents, or inter- 
action with its environment (the "thermal bath"), etc. Following Gardiner 
[2] we shall give a semi-historical outline of how a phenomenological theory 
of fluctuating phenomena arose and what its essential points are. 

The experience of careful measurements in science normally gives us data 
like that of Fig. [U representing the time evolution of a certain variable 
X. Here a quite well defined deterministic trend is evident, which is re- 
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Figure 1: Schematic time evolution of a variable X with a well defined 
deterministic motion plus fluctuations around it. 
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producible, unlike the fluctuations around this motion, which are not. This 
evolution could represent, for instance, the growth of the (normalised) num- 
ber of molecules of a substance X formed by a chemical reaction of the form 
A ^ X, or the process of charge of a capacitor in a electrical circuit, etc. 

1.1 Brownian motion 

The observation that, when suspended in water, small pollen grains are found 
to be in a very animated and irregular state of motion, was first systemati- 
cally investigated by Robert Brown in 1827, and the observed phenomenon 
took the name of Brownian motion. This motion is illustrated in Fig. [2J 
Being a botanist, he of course tested whether this motion was in some way 
a manifestation of life. By showing that the motion was present in any sus- 
pension of fine particles — glass, mineral, etc. — he ruled out any specifically 
organic origin of this motion. 

1.1.1 Einstein's explanation (1905) 

A satisfactory explanation of Brownian motion did not come until 1905, when 
Einstein published an article entitled Concerning the motion, as required 
by the molecular-kinetic theory of heat, of particles suspended in liquids at 
rest. The same explanation was independently developed by Smoluchowski 




Figure 2: Motion of a particle undergoing Brownian motion. 
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in 1906, who was responsible for much of the later systematic development 
of the theory. To simplify the presentation, we restrict the derivation to a 
one-dimensional system. 

There were two major points in Einstein's solution of the problem of 
Brownian motion: 

• The motion is caused by the exceedingly frequent impacts on the pollen 
grain of the incessantly moving molecules of liquid in which it is sus- 
pended. 

• The motion of these molecules is so complicated that its effect on the 
pollen grain can only be described probabilistically in term of exceed- 
ingly frequent statistically independent impacts. 

Einstein development of these ideas contains all the basic concepts which 
make up the subject matter of these notes. His reasoning proceeds as follows: 
"It must clearly be assumed that each individual particle executes a motion 
which is independent of the motions of all other particles: it will also be 
considered that the movements of one and the same particle in different time 
intervals are independent processes, as long as these time intervals are not 
chosen too small." 

"We introduce a time interval r into consideration, which is very small 
compared to the observable time intervals, but nevertheless so large that in 
two successive time intervals r, the motions executed by the particle can be 
thought of as events which are independent of each other." 

"Now let there be a total of n particles suspended in a liquid. In a time 
interval r, the X-coordinates of the individual particles will increase by an 
amount A, where for each particle A has a different (positive or negative) 
value. There will be a certain frequency law for A; the number dn of the 
particles which experience a shift between A and A + dA will be expressible 
by an equation of the form: dn = n0(A)dA, where 0(A)dA = 1, and 
is only different from zero for very small values of A, and satisfies the 
condition 0(-A) = 0(A)." 

"We now investigate how the diffusion coefficient depends on 0. We shall 
restrict ourselves to the case where the number of particles per unit volume 
depends only on x and t." 

"Let f(x,t) be the number of particles per unit volume. We compute 
the distribution of particles at the time t + r from the distribution at time 
t. From the definition of the function 0(A), it is easy to find the number of 
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particles which at time t + r are found between two planes perpendicular to 
the a>axis and passing through points x and x + dx. One obtains: 




(1.1) 



But since r is very small, we can set 



df 

Furthermore, we expand f(x + A, t) in powers of A: 

* s s ,df(x,t) A 2 d 2 f(x,t) 

We can use this series under the integral, because only small values of A 
contribute to this equation. We obtain 
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Because </>(— A) = 0(A), the second, fourth, etc. terms on the right-hand 
side vanish, while out of the 1st, 3rd, 5th, etc., terms, each one is very small 
compared with the previous. We obtain from this equation, by taking into 
consideration 

/oo 
0(A)dA = 1 . 
-oo 

and setting 

1 f°° A 2 

-0(A)dA = D , (1.3) 



rj-oo 2 



and keeping only the 1st and 3rd terms of the right hand side, 



dt dx 2 



(1.4) 



This is already known as the differential equation of diffusion and it can be 
seen that D is the diffusion coefficient." 
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Figure 3: Time evolution of the non-equilibrium probability distribution 
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"The problem, which correspond to the problem of diffusion from a sin- 
gle point (neglecting the interaction between the diffusing particles), is now 
completely determined mathematically: its solution is 



f(x,t) 



-x 2 /4Dt 



;i.s) 



This is the solution, with the initial condition of all the Brownian particles 
initially at x = 0; this distribution is shown in Fig. [3] 



1 We can get the solution (|1.5[) by using the method of the integral transform to solve 
partial differential equations. Introducing the space Fourier transform of f(x,t) and its 
inverse, 

F(k, t)= J Ax e~ ikx f(x, t) , f{x, t) = ±-j&k e ikx F(k, t) , 
the diffusion equation (|1.4p transforms into the simple form 



dt 



F(k,t) = F{k,0)e 



-Dk J t 



For the initial condition f(x,t = 0) = S(x), the above Fourier transform gives F(k,t 
0) = 1. Then, taking the inverse transform of the solution in fc-space, we finally have 



f(x,t) = ±- /.!/,-< ■ LU - li 

Z7T 



„-x 2 /4Dt r p -x 2 /4Dt 

1 dke- Dt( - k ~ ix/2Dt) 



2n 



where in the second step we have completed the square in the argument of the exponential 
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Einstein ends with: "We now calculate, with the help of this equation, 
the displacement X x in the direction of the X-axis that a particle experiences 
on the average or, more exactly, the square root of the arithmetic mean of 
the square of the displacements in the direction of the X-axis; it is 



Einstein derivation contains very many of the major concepts which since 
then have been developed more and more generally and rigorously over the 
years, and which will be the subject matter of these notes. For example: 

(i) The Chapman-Kolgomorov equation occurs as Eq. (II. II) . It states that 
the probability of the particle being at point x at time t + r is given by 
the sum of the probabilities of all possible "pushes" A from positions 
x + A, multiplied by the probability of being at x + A at time t. 
This assumption is based on the independence of the push A of any 
previous history of the motion; it is only necessary to know the initial 
position of the particle at time t — not at any previous time. This is the 
Markov postulate and the Chapman-Kolmogorov equation, of which 
Eq. fll.il) is a special form, is the central dynamical equation to all 
Markov processes. These will be studied in Sec. |3j 

(ii) The Kramer s-Moyal expansion. This is the expansion used [Eq. (|1.2)) ] 
to go from Eq. (11.11) (the Chapman-Kolmogorov equation) to the dif- 
fusion equation (11.41) . 

(iii) The Fokker-Planck equation. The mentioned diffusion equation (jl.4p . 
is a special case of a Fokker-Planck equation. This equation governs 
an important class of Markov processes, in which the system has a 
continuous sample path. We shall consider points (ii) and (iii) in detail 
in Sec. SJ 

1.1.2 Langevin's approach (1908) 

Some time after Einstein's work, Langevin presented a new method which 
was quite different from the former and, according to him, "infiniment plus 
simple" . His reasoning was as follows. 

—Dk 2 t + ikx = —Dt(k — \x/2Dt) 2 — x 2 /4Dt, and in the final step we have used the 
Gaussian integral Jdk e - a ( k ~ b ) — yV/a, which also holds for complex b. 




(1.6) 
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From statistical mechanics, it was known that the mean kinetic energy of 
the Brownian particles should, in equilibrium, reach the value 



{\mv 2 ) = \k B T . (1.7) 
Acting on the particle, of mass m, there should be two forces: 

(i) a viscous force: assuming that this is given by the same formula as 
in macroscopic hydrodynamics, this is — rwydx/dt, with 7727 = 6irfia, 
being /i the viscosity and a the diameter of the particle. 

(ii) a fluctuating force £(£), which represents the incessant impacts of the 
molecules of the liquid on the Brownian particle. All what we know 
about it is that is indifferently positive and negative and that its mag- 
nitude is such that maintains the agitation of the particle, which the 
viscous resistance would stop without it. 

Thus, the equation of motion for the position of the particle is given by 
Newton's law as 

(\ 2 T At 

fl. 




Multiplying by x, this can be written 

md 2 (i 2 ) 2 772,7 d(a; 2 ) 
2 dt 2 2 dt s 

If we consider a large number of identical particles, average this equation 
written for each one of them, and use the equipartition result (II. 7p for (mv 2 ), 
we get and equation for (x 2 ) 

m d 2 (x 2 ) 7777 d (x 2 ) 
~2^t 2 - + —— = kBT - 

The term (£ x) has been set to zero because (to quote Langevin) "of the irreg- 
ularity of the quantity £(i)" . One then finds the solution (C is an integration 
constant) 

d (r 2 ) 

= 2k B T/ mi + Ce~^ . 

dt 
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Langevin estimated that the decaying exponential approaches zero with a 
time constant of the order of 10~ 8 s, so that d (x 2 ) /dt enters rapidly a con- 
stant regime d (x 2 ) /dt = Ik^T/m^ Therefore, one further integration (in 
this asymptotic regime) leads to 



which corresponds to Einstein result (11.61) . provided we identify the diffusion 
coefficient as 



It can be seen that Einstein's condition of the independence of the displace- 
ments A at different times, is equivalent to Langevin's assumption about the 
vanishing of (£ x). Langevin's derivation is more general, since it also yields 
the short time dynamics (by a trivial integration of the neglected Ce~ 7 *), 
while it is not clear where in Einstein's approach this term is lost. 

Langevin's equation was the first example of a stochastic differential equa- 
tion — a differential equation with a random term £(t) and hence whose so- 
lution is, in some sense, a random functional Each solution of the Langevin 
equation represents a different random trajectory and, using only rather 
simple properties of the fluctuating force £(£), measurable results can be 
derived. Figure H] shows the trajectory of a Brownian particle in two dimen- 
sions obtained by numerical integration of the Langevin equation (we shall 
also study numerical integration of stochastic differential equations). It is 
seen the growth with t of the area covered by the particle, which corresponds 
to the increase of (x 2 ) — (x^) in the one-dimensional case discussed above. 

The theory and experiments on Brownian motion during the first two 
decades of the XX century, constituted the most important indirect evidence 
of the existence of atoms and molecules (which were unobservable at that 
time). This was a strong support for the atomic and molecular theories of 
matter, which until the beginning of the century still had strong opposition 
by the so-called energeticits. The experimental verification of the theory of 
Brownian motion awarded the 1926 Nobel price to Svedberg and Perrin. [f| 

2 The rigorous mathematical foundation of the theory of stochastic differential equa- 
tions was not available until the work of Ito some 40 years after Langevin's paper. 

3 Astonishingly enough, the physical basis of the phenomenon was already described 
in the 1st century B.C.E. by Lucretius in De Rerum Natura (II, 112-141), a didactical 
poem which constitutes the most complete account of ancient atomism and Epicureanism. 



(x 2 ) - (x 2 ) = 2{k B T/ mi )t , 



D = k-BT/rwy . 
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Figure 4: Trajectory of a simulated Brownian particle projected into the 
x-y plane, with D = 0.16 /zm 2 /s. The x and y axes are marked in microns. 
It starts from the origin (x,y) = (0,0) at t — 0, and the pictures show the 
trajectory after 1 sec, 3 sec and 10 sec. 
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The picture of a Brownian particle immersed in a fluid is typical of a 
variety of problems, even when there are no real particles. For instance, it 
is the case if there is only a certain (slow or heavy) degree of freedom that 
interacts, in a more or less irregular or random way, with many other (fast 
or light) degrees of freedom, which play the role of the bath. Thus, the 
general concept of fluctuations describable by Fokker-Planck and Langevin 
equations has developed very extensively in a very wide range of situations. 
A great advantage is the necessity of only a few parameters; in the example 
of the Brownian particle, essentially the coefficients of the derivatives in the 
Kramers-Moyal expansion (allowing in general the coefficients a x and t 
dependence) 



It is rare to find a problem (mechanical oscillators, fluctuations in electri- 
cal circuits, chemical reactions, dynamics of dipoles and spins, escape over 
metastable barriers, etc.) with cannot be specified, in at least some degree 
of approximation, by the corresponding Fokker-Planck equation, or equiv- 
alently, by augmenting a deterministic differential equation with some fluc- 
tuating force or field, like in Langevin's approach. In the following sections 
we shall describe the methods developed for a systematic and more rigorous 
study of these equations. 



When observing dust particles dancing in a sunbeam, Lucretius conjectured that the 
particles are in such irregular motion since they are being continuously battered by the 
invisible blows of restless atoms. Although we now know that such dust particles' motion 
is caused by air currents, he illustrated the right physics but only with a wrong example. 
Lucretius also extracted the right consequences from the "observed" phenomenon, as one 
that shows macroscopically the effects of the "invisible atoms" and hence an indication of 
their existence. 




(1.10) 
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2 Stochastic variables 



2.1 Single variable case 

A stochastic or random variable is a quantity X, denned by a set of possible 
values {x} (the "range", "sample space", or "phase space"), and a probability 
distribution on this set, Px-(^)0 The range can be discrete or continuous, 
and the probability distribution is a non-negative function, P x {x) > 0, with 
P x (x)dx the probability that X £ (x, x + dx). The probability distribution 
is normalised in the sense 

JdxP x (x) = l, 

where the integral extends over the whole range of X. 

In a discrete range, {x n }, the probability distribution consists of a number 
of delta-type contributions, P x { x ) — J2 n Pn8(x — x n ) and the above normal- 
isation condition reduces to ^2 n p n = 1- For instance, consider the usual 
example of casting a die: the range is {x n } = {1,2,3,4,5,6} and p n = 1/6 
for each x n (in a honest die). Thus, by allowing 5-function singularities in 
the probability distribution, one may formally treat the discrete case by the 
same expressions as those for the continuous case. 

2.1.1 Averages and moments 

The average of a function f(X) defined on the range of the stochastic variable 
X, with respect to the probability distribution of this variable, is defined as 

(f{X)) = Jdxf{x)P x {x) . 

The moments of the stochastic variable, /i m , correspond to the special cases 
f(X)=X m ,i.e$ 

/x m = (X m ) = Jdxx m P x (x) , m = l,2,.... (2.1) 

4 It is advisable to use different notations for the stochastic variable, X, and for the 
corresponding variable in the probability distribution function, x. However, one relaxes 
this convention when no confusion is possible. Similarly, the subscript X is here and there 
dropped from the probability distribution. 

5 This definition can formally be extended to m — 0, with /i = 1, which expresses the 
normalisation of Px (x) . 
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2.1.2 Characteristic function 

This useful quantity is defined by the average of exp(ikX), namely 

G x (k) = (exp(iJfeX)) = Jdx exp(ikx)P x (x) . (2.2) 

This is merely the Fourier transform of Px(x), and can naturally be solved 
for the probability distribution 



Px(x) 



1 [°° 

— / dk exp(—ikx)Gx(k) . 
2vr i_oo 



The function Gx{k) provides an alternative complete characterisation of the 
probability distribution. 

By expanding the exponential in the integrand of Eq. (12.21) and inter- 
changing the order of the resulting series and the integral, one gets 

Gx{k) = / dzx-P^x) = V ■ (2.3) 

z — ' ml J ml 

m=0 m=0 

Therefore, one finds that Gx{k) is the moment generating function, in the 
sense that 

»™ = ("0 m -^Gx(k) . 



dk r 



(2.4) 

fc=0 



2.1.3 Cumulants 



The cumulants, K m , are defined as the coefficients of the expansion of the 
cumulant function hiGx(k) in powers of ik, that is, 

°^ (i£;) m 

\nG x {k) = V 7-n m . 

/ — ' ml 

m=l 

Note that, owing to Px{x) is normalised, the m = term vanishes and the 
above series begins at m = 1. The explicit relations between the first four 
cumulants and the corresponding moments are 



K-i = /ii 

«2 = fJ<2 - lA 



^3 = Ai3 - 3/i 2 /ii + 2/if 
k 4 = /x 4 - 4/i 3/ ui - 3/x| + 12/i 2 /i? - 6/4 ■ 
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Thus, the first cumulant is coincident with the first moment (mean) of the 
stochastic variable: K\ = (X); the second cumulant K2, also called the 
variance and written a 2 , is related to the first and second moments via 
a 2 = k 2 = We finally mention that there exists a general 

expression for K m in terms of the determinant ofamxm matrix constructed 
with the moments {/Xj | 



\m— 1 



1, 


. . ,m} 


(see, 


e.g., H p. 18]) 
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(2.6) 



where the Q) are binomial coefficients. 



2.2 Multivariate probability distributions 

All the above definitions, corresponding to one variable, are readily extended 
to higher- dimensional cases. Consider the n-dimensional vector of stochastic 
variables X = (Xi, . . . , X n ), with a probability distribution P n (xi, . . . , x n ). 
This distribution is also referred to as the joint probability distribution and 

Pni^Xli • • • j X n )dXi • • • dx n , 

is the probability that X±, . . . , X n have certain values between (x±,xi + 
dxi),. . . , {x n ^ x n -\- dx n ). 



Partial distributions. One can also consider the probability distribution 
for some of the variables. This can be done in various ways: 

1. Take a subset of s < n variables X±, . . . , X s . The probability that they 
have certain values in (xi, X\ + dxi), . . . , (x s , x s + dx s ), regardless of the 
values of the remaining variables X s+ i, . . . , X n , is 

Psip^Xt • • • 1 Xs) / d-Xg+i • • • dx n P n (xi, . . . , X s , . . . , x n ) , 



6 Quantities related to the third- and fourth-order cumulants have also their own names: 
skewness, k 3 /k^ 2 , and kurtosis, K4/ ' n\. 
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which is called the marginal distribution for the subset Xi, . . . ,X S . 
Note that from the normalisation of the joint probability it immediately 
follows the normalisation of the marginal probability. 

2. Alternatively, one may attribute fixed values to X s+1 , . . . , X n , and con- 
sider the joint probability of the remaining variables Xi, . . . , X s . This 
is called the conditional probability, and it is denoted by 

Ps\n—s{X\ i ■ ■ ■ i X s \ x s +i, . . . , Xn) . 

This distribution is constructed in such a way that the total joint proba- 
bility P n (xi, ■ ■ ■ , x n ) is equal to the marginal probability for X s+ ±, . . . , X n 
to have the values x s+ i, . . . ,x n , times the conditional probability that, 
this being so, the remaining variables X±, . . . , X s have the values (xi, . . . , x s ). 
This is Bayes' rule, and can be considered as the definition of the con- 
ditional probability: 

Pn(x\ ■>•••■> X n ) Pn— s(Xs+li ■ ■ ■ i Xn)P s \ n — s {x\ , • • • , X s \x s +\ , . . . , X n ) . 

Note that from the normalisation of the joint and marginal probabilities 
it follows the normalisation of the conditional probability. 



Characteristic function: moments and cumulants. For multivariate 
probability distributions, the moments are defined by 

(Xr-.-X™") = Jdx 1 ---dx n xT 1 ---x^P(x 1 ,...,x n ) , 

while the characteristic (moment generating) function depends on n auxiliary 
variables k = (k± : . . . , k n ): 



G(k) = (exp[i(k 1 X 1 + --- + k n X n )}) 

l) mi ■ ■ ■ (ifcn 

• ■ ■ m n \ 





Similarly, the cumulants of the multivariate distribution, indicated by double 
brackets, are defined in terms of the coefficients of the expansion of In G as 

lnG(k)=£' ^ ^ ((Xr-.-X^)) , 

Ul\. • ■ ■ IIL n . 

where the prime indicates the absence of the term with all the m 8 simulta- 
neously vanishing (by the normalisation of P n ). 
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Covariance matrix. The second-order moments may be combined into a 
n by n matrix (XiXj). More relevant is, however, the covariance matrix, 
defined by the second-order cumulants 



Each diagonal element is called the variance of the corresponding variable, 
while the off-diagonal elements are referred to as the covariance of the cor- 
responding pair of variables^ 

Statistical independence. A relevant concept for multivariate distribu- 
tions is that of statistical independence. One says that, e.g., two stochastic 
variables X\ and X 2 are statistically independent of each other if their joint 
probability distribution factorises: 



The statistical independence of X\ and X 2 is also expressed by any one of 
the following three equivalent criteria: 

1. All moments factorise: (X^X™ 2 ) = (X™ 1 ) (X™ 2 ) . 

2. The characteristic function factorises: Gi 1 x 2 (^i ) k 2 ) = G x^k^G x 2 ik 2 ) • 

3. The cumulants ((X^X™ 2 )) vanish when both mi and m 2 are ^ 0. 

Finally, two variables are called uncorrelated when its covariance, ((XiX 2 )), 
is zero, which is a condition weaker than that of statistical independence. 

2.3 The Gaussian distribution 

This important distribution is defined as 



It is easily seen that \L\ is indeed the average and a 1 the variance, which 
justifies the notation. The corresponding characteristic function is 






(2.9) 



7 Notc that ((XiXj)) is, by construction, a symmetrical matrix. 
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as can be seen from the definition (12.21) . by completing the square in the 
argument of the total exponential ikx — (x — fii) 2 /2a 2 and using the Gaussian 
integral J dk e~ a ( k ~ b ^ = \/ir/a for complex b as in the footnote in p. [7J 
Note that the logarithm of this characteristic function comprises terms up 
to quadratic in k only. Therefore, all the cumulants after the second one 
vanish identically, which is a property that indeed characterises the Gaussian 
distribution. 

For completeness, we finally write the Gaussian distribution for n vari- 
ables X = (Xl, . . . , X n ) and the corresponding characteristic function 

-i(x -x ) • A ■ (x - x ) 
G(k) = exp (ix • k - ik ■ i- 1 • k) . 

The averages and covariances are given by (X) = x and ((JQJfj)) = (A^ 1 )^. 



PU) 



aei i\ 
(2ir) n 



exp 



2.4 Transformation of variables 

For a given stochastic variable X, every related quantity Y = f(X) is again 
a stochastic variable. The probability that Y has a value between y and 
V + Ay is 

P Y {y)Ay= JdxP x {x), 

y<f(x)<y+Ay 

where the integral extends over all intervals of the range of X where the 
inequality is obeyed. Note that one can equivalently define Py(y) as^l 




(2.11) 



8 Note also that from Eq. (|2.11|1 . one can formally write the probability distribution 
for Y as the following average [with respect to Px (x) and taking y as a parameter] 



Py(y) = (S[y-f(X)}) 



(2.10) 
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From this expression one can calculate the characteristic function of Y: 



G Y (k) 



Eq. 




dy exp(iky)P Y (y) 

dxPx(x) J dy exjp(iky)6[y - f(x)) 

dxP.Wexp^W], 



which can finally be written as 



G Y (k) = <exp[ifc/(X)]> . 



(2.12) 



As the simplest example consider the linear transformation Y = aX. The 
above equation then yields Gy{k) = (exp(ifcaX)), whence 



2.5 Addition of stochastic variables 

The above equations for the transformation of variables remain valid when 
X stands for a stochastic variable with n components and Y for one with 
s components, where s may or may not be equal to n. For example, let us 
consider the case of the addition of two stochastic variables Y = f(Xi, X2) = 
X\ + X2, where s — 1 and n — 2. Then, from Eq. (12.111) one first gets 



Properties of the sum of stochastic variables. One easily deduces the 
following three rules concerning the addition of stochastic variables: 



G Y (k) = G x {ak) 



(Y = aX) . 



(2.13) 




dxi dx 2 5[y - (xx + x 2 )]Px 1 x 2 (x 1 ,x 2 ) 



(2.14) 
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1. Regardless of whether X\ and X 2 are independent or not, one 

(Y) = (X 1 ) + (X 2 ) . (2.15) 

2. If X\ and X 2 are uncorrelated, ((X1X2)) = 0, a similar relation holds 
for the variances 10 ! 



«^» = «*?» + «*f» ■ (2-W) 
3. The characteristic function of Y = X\ + X 2 is^l 

Gy(A;) = G Xl x 2 (k, k) . (2.17) 

On the other hand, if Xi and X 2 are independent, Eq. (12.141) and the factor- 
ization of P Xl x 2 an d ^^1X2 yields 

P Y (y) = [ dx 1 P Xl (x 1 )P X2 (y-x 1 ) , Gy(A;) = G Xl {k)G X2 {k) . (2.18) 



9 Proof of Eq. ([2715]) : 

(F) = JdyyP Y (y) = Jdx x J dx 2 P Xl x 2 {xi, x 2 ) J dy yS[y - {x x + x 2 )] 
= Jdxx Jdx 2 Px 1 x 2 (xi,x 2 )(x 1 + x 2 ) = (Jfi) + (X 2 ) ■ Q.E.D. 

10 Proof of Eq. ([2715]) : 

(F 2 ) = /" dyy 2 /V(y) = /"dari / dx 2 Px^ (an, a; 2 )(aji +x 2 ) 2 = (JC 1 2 ) + (X 2 2 )+2 {X X X 2 ) 
Therefore 

((y 2 )) = (y 2 ) - (Yf = (xf) (x,) 2 + (xi) (x 2 f +2 {(x x x 2 ) (x 1 ) (x 2 )) 

((x ? )) ((xi)) 

from which the statement follows for uncorrelated variables. Q.E.D. 

11 Proof of Eq. (|2~I7|) : 

Gy(fc) ^i 111 (exp[ifc/(Xi,Jf2)]) - (exp[ife(Xi + X 2 )]) G Xl x 2 (k,k) . Q.E.D. 
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Thus, the probability distribution of the sum of two independent variables 
is the convolution of their individual probability distributions. Correspond- 
ingly, the characteristic function of the sum [which is the Fourier transform 
of the probability distribution; see Eq. (12.21) ] is the product of the individual 
characteristic functions. 



2.6 Central limit theorem 

As a particular case of transformation of variables, one can also consider the 
sum of an arbitrary number of stochastic variables. Let X 1} . . . ,X n be a 
set of n independent stochastic variables, each having the same probability 
distribution Px{x) with zero average and (finite) variance a 2 . Then, from 
Eqs. ( 12. 15)) and (12.161) it follows that their sum Y — X\ -\ — • + X n has zero 
average and variance no~ 2 , which grows linearly with n. On the other hand, 
the distribution of the arithmetic mean of the variables, (Xi + ■ • • + X n )/n, 
becomes narrower with increasing n (variance a 2 /n). It is therefore more 
convenient to define a suitable scaled sum 

Xl + ..- + x n 

~ T ' 

which has variance a 2 for all n. 

The central limit theorem states that, even when Px{x) is not Gaus- 
sian, the probability distribution of the Z so-defined tends, as n —>■ oo, to 
a Gaussian distribution with zero mean and variance a 2 . This remarkable 
result is responsible for the important role of the Gaussian distribution in all 
fields in which statistics are used and, in particular, in the equilibrium and 
non-equilibrium statistical physics. 



Proof of the central limit theorem. We begin by expanding the char- 
acteristic function of an arbitrary Px{x) with zero mean as [cf. Eq. (12. 3p 

G x {k) = Jdx exp(ikx)P x (x) = 1 - \a 2 k 2 + ■■■ . (2.19) 

The factorization of the characteristic function of the sum Y = Xi + • • • + X n 
of statistically independent variables [Eq. (12.181) ]. yields 

n 

G Y (k) = l[GxAk) = [Gx(k)} n , 



i=l 
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where the last equality follows from the equivalent statistical properties of 
the different variables Xj. Next, on accounting for Z = Y/ y/n, and using the 
result (I2.13P with a = 1 / y/n, one has 

°^ - ^ ty - h (£)]" * l 1 - ™ exp (4<r2fc2) ' (2 - 2o) 

where we have used the definition of the exponential e x = lim (1 + x/n) n . 

n— >oo 

Finally, on comparing the above result with Eqs. (12.81) . one gets 



P z (z) = ^L=exp(-^—) . Q.E.D. 



. 2 



Remarks on the validity of the central limit theorem. The derivation 
of the central limit theorem can be done under more general conditions. 
For instance, it is not necessary that all the cumulants (moments) exist. 
However, it is necessary that the moments up to at least second order exist 
[or, equivalently, Gx{k) being twice differentiable at k = 0; see Eq. (12.41) . 
The necessity of this condition is illustrated by the counter-example provided 
by the Lorentz-Cauchy distribution: 

1 7 

P(x) = = , (— oo < x < oo) . 

7T X 1 + Y 

It can be shown that, if a set of n independent variables Xj have Lorentz- 
Cauchy distributions, their sum also has a Lorentz-Cauchy distribution (see 
footnote below). However, for this distribution the conditions for the cen- 
tral limit theorem to hold are not met, since the integral (12.11) defining the 
moments yU m , does not converge even for m = ljjfl 

Finally, although the condition of independence of the variables is impor- 
tant, it can be relaxed to incorporate a sufficiently weak statistical depen- 
dence. 



12 This can also be demonstrated by calculating the corresponding characteristic func- 
tion. To do so, one can use f_ °°dx e lax /(l + x 2 ) — 7re~! a ', which is obtained by comput- 
ing the residues of the integrand in the upper (lower) half of the complex plane for a > 
(a < 0). Thus, one gets 

G(k) = exp(- 7 |fc|) , 

which, owing to the presence of the modulus of k, is not differentiable at k — 0. Q.E.D. 

We remark in passing that, from Gxt(k) = exp(— 7i|fc|) and the second Eq. (|2. 18|) . it 
follows that the distribution of the sum of independent Lorentz-Cauchy variables has a 
Lorentz-Cauchy distribution (with Gy{k) = expj— Q^. 7i)|fc|])- 
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2.7 Exercise: marginal and conditional probabilities 
and moments of a bivariate Gaussian distribution 



To illustrate the definitions given for multivariate distributions, let us com- 
pute them for a simple two-variable Gaussian distribution 



P 2 {xi,x 2 ) 



1 - A 2 
^W ex P 



(27T(T 2 



2a 2 



x. 



2A x\Xi + x 2 ,) 



^2.23) 



where A is a parameter — 1 < A < 1, to ensure that the quadratic form in 
the exponent is definite positive (the equivalent condition to assume a 2 to be 
positive in the one- variable Gaussian distribution (12. 8p . The normalisation 
factor can be seen to take this value by direct integration, or by comparing 
our distribution with the multidimensional Gaussian distribution (Sec. 12. 3j) : 

here A = — ^ jj^ ^ j so that det A = (1 — A 2 )/cr 4 . Finally, if one wishes 

to fix ideas one can interpret P 2 (xi, x 2 ) as the Boltzmann distribution of two 
harmonic oscillators coupled by a potential term oc XxiX 2 - 

Let us first rewrite the distribution in a form that will facilitate to do the 
integrals by completing once more the square — 2A X\X 2 + x 2 = (x 2 — Axi) 2 — 
A 2 x 2 



P 2 (xi, x 2 ) = C exp 







2a 2 % \ 


exp 



2a 2 



(x 2 - Axi) 2 



^2.24) 



and C = \/ (1 — A 2 ) / (27rcr 2 ) 2 is the normalisation constant. We can now 
compute the marginal probability of the individual variables (for one of them 
since they are equivalent), defined by P\{x\) = J dx 2 P 2 (xi,x 2 ) 

p^xx) = Ce~^ x * [dx 2 e- {X2 - Xxi)2/2a2 . 



n/a a=l/ 2 a 2 



Therefore, recalling the form of C, we merely have 



Pi{xi) 



/ 2tto 2 



exp 



xA 



with o\ = a 2 /(l - A 2 ) . (2.25) 



We see that the marginal distribution depends on A, which results in a 
modified variance. To see that a 2 is indeed the variance ((x 2 )) = (x 2 ) — (xi) 2 , 
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note that {x™) can be obtained from the marginal distribution only (this is 
a general result) 

{xf)= dxt dx 2 xfB 2 (x 1: x 2 ) = I dx x x™ I dx 2 B 2 {x u x 2 ) = / dx x xTP\{x\) 



Pi(xi) 



Then inspecting the marginal distribution obtained [Eq. f)2.25p ] we get that 
the first moments vanish and the variances are indeed equal to a 



2. 




(2.26) 



To complete the calculation of the moments up to second order we need 
the covariance of x\ and x 2 . {{x\x 2 )) = (x\x 2 ) — (xi) (x 2 ) which reduces 
to calculate {x\X 2 ). This can be obtained using the form (I2.24p for the 
distribution 



(XiX 2 ) 



J dx 1 J dx 2 XiX 2 P 2 (x 1 ,x 2 ) 



C / dxiXi exp 



1 - A 



2(7 



2 1 



exp 




• v /27r CT 2 /(l-A 2 )(T 2 /(l-> 2 ) 



since C = yJJY— A 2 )/(27rcr 2 ) 2 . Its is convenient to compute the normalised 
covariance (xix 2 ) / a/ (xf) which is merely given by A. Therefore the 
parameter A in the distribution f!2.23j) is a measure of how much correlated 
the variables x\ and x 2 are. Actually in the limit A —>■ the variables are 
not correlated at all and the distribution factorises. In the opposite limit 
A — > 1 the variables are maximally correlated, (xix 2 ) / a/ (xf) (x 2 ) = 1. The 
distribution is actually a function of (x± — x 2 ), so it is favoured that x\ and 
x 2 take similar values (see Fig. [5]) 



2|A=0 



-xl/2a 2 _ 



-x\l2a 2 



V2 



V2 



2 A=l 



-( Xl - X2 )l/2a2 



. (2.27) 



7T& 
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We can now interpret the increase of the variances with A: the correlation 
between the variables allows them to take arbitrarily large values, with the 
only restriction of their difference being small (Fig. [5]). 

To conclude we can compute the conditional probability by using Bayes 
rule Pi|i(xi|x 2 ) = P 2 (xi,x 2 )/Pi(x 2 ) and Eqs. (12351) and ([2^25]) 



Pi\i(x 1 \x 2 ) 



l-A 2 



exp [—;^2 {%\ — 2\xiX 2 + x 2 )] 



l-A 2 



V2 



: CXp 



2^ eX P( 
1 



l-A 2 J2 



2a 2 







2(T 2 



(xl-2Xx l x 2 + [y-(y-x 2 )]x 2 2 ) 



exp(-0.5*(x"2-2.'0.*x*y+y"2)) exp(-0.5-(x"2-2."0.6'x-y+y"2)) 

0.9 0.9 




exp(-0.5-(x"2-2.-D.9'x-y+y"2)| - 




exp(-0.5 - (x"2-2.*1 .'x'y+y"2)) - 




Figure 5: Gaussian distribution (I2.23I) for A = 0, 0.6, 0.9 and 1 (non- 
normalised) . 
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and hence (recall that here x 2 is a parameter; the known output of X 2 ) 



Pi\i(xi\x 2 ) 



V2 



: CXp 



7T(J Z 



2a 2 



(2.28) 



Then, at A = (no correlation) the values taken by x\ are independent of 
the output of x 2 while for A — > 1 they are centered around those taken by 
x 2 , and hence strongly conditioned by them. 
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3 Stochastic processes and Markov processes 



Once a stochastic variable X has been defined, an infinity of other stochastic 
variables derive from it, namely, all the quantities Y defined as functions of 
X by some mapping. These quantities Y may be any kind of mathematical 
object; in particular, also functions of an auxiliary variable t: 

Y x (t)=f(X,t), 

where t could be the time or some other parameter. Such a quantity Yx(t) 
is called a stochastic process. On inserting for X one of its possible values x, 
one gets an ordinary function of t, Y x (t) = f(x,t), called a sample function 
or realisation of the process. In physical language, one regards the stochastic 
process as the "ensemble" of these sample functions 

It is easy to form averages on the basis of the underlying probability 
distribution Px(x). For instance, one can take n values t\,...,t n , and form 
the nth moment 

(Y(t 1 )---Y(t n )) = jdxY x {t x )---Y x {t n )P x {x) . (3.1) 

Of special interest is the auto- correlation function 

K(t u t 2 ) = ((Y( tl )Y(t 2 )}} = (Y( tl )Y(t 2 )) - (Y( tl )) (Y(t 2 )) 

= ([Y^-iY^MY^-iYit,))]), 

which, for t\ — t 2 — t, reduces to the time-dependent variance ((Y(t) 2 )) = 
a{tf. 

A stochastic process is called stationary when the moments are not af- 
fected by a shift in time, that is, when 

(y(ti + r) • ■ • Y(t n +t)) = (Y(tx) ■ ■ ■ Y(t n )} . (3.2) 

In particular, (Y(t)) is then independent of the time, and the auto-correlation 
function K(t\,t 2 ) only depends on the time difference \ti — t 2 \. Often there 
exist a constant r c such that n(t\,t 2 ) ~ for \t\ — t 2 \ > r c ; one then calls r c 
the auto- correlation time of the stationary stochastic process. 

13 As regards the terminology, one also refers to a stochastic time-dependent quantity as 
a noise term. This name originates from the early days of radio, where the great number 
of highly irregular electrical signals occurring either in the atmosphere, the receiver, or 
the radio transmitter, certainly sounded like noise on a radio. 
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If the stochastic quantity consist of several components Yi(t), the auto- 
correlation function is replaced by the correlation matrix 

K^h) = ((Y^Yjih))) . 

The diagonal elements are the auto- correlations and the off-diagonal ele- 
ments are the cross-correlations. Finally, in case of a zero-average stationary 
stochastic process, this equation reduces to 

K t3 {r) = (Y(t)Y,(t + r)) = (Yj(O)Yj-(T)) . 

3.1 The hierarchy of distribution functions 

A stochastic process Yx(t), defined from a stochastic variable X in the way 
described above, leads to a hierarchy of probability distributions. For in- 
stance, the probability distribution for Yx(t) to take the value y at time t is 
[cf. Eq. fl2TU)] 

Pi{y,t) = Jdx5[y-YJ^\P x (x) . 

Similarly, the joint probability distribution that Y has the value y% at t±, and 
also the value 2/2 at and so on up to y n at t n , is 

P n (Vu W, O = |dz%, - YAh)] ■ ■ ■ %„ - YM\P X (X) . (3.3) 

In this way an infinite hierarchy of probability distributions P n , n = 1, 2, . . ., 
is defined. They allow one the computation of all the averages already intro- 
duced, e.g.0 

(Y(ti) ■ ■ - Y(t n )) = dyi---dy n yx---y n P n (y 1 ,ti;...;y n ,t n ). (3.4) 



14 This result is demonstrated by introducing the definition (|3.3|) in the right-hand side 
of Eq. (P0j> : 

Jdyi ■ ■ ■ dy n y 1 ■ ■ ■ y n P n (vuh; ■■■] Vn, t n ) 

dyi • ■■dy n yx ■■■y n / dx<5[yi - Y x (tx)] ■ ■ ■ S[y n - Y x (t n )]P x {x) 



Eq. d37TJ 



(Y(t 1 )---Y(t n )) . Q.E.D. 
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We note in passing that, by means of this result, the definition (13. 2p of 
stationary processes, can be restated in terms of the dependence of the P n 
on the time differences alone, namely 

P n (yi, h + T;...;y n ,t n + T) = P n (yi, h;...;y n , t n ) . 

Consequently, a necessary (but not sufficient) condition for the stochastic 
process being stationary is that P\{yi) does not depend on the time. 

Although the right-hand side of Eq. (I3.3P also has a meaning when some 
of the times are equal, one regards the P n to be defined only when all times 
are different. The hierarchy of probability distributions P n then obeys the 
following consistency conditions: 

1. P n > o ; 

2. P n is invariant under permutations of two pairs (?/;,£;) and (yj,tj) ; 

3. /d V „P„( !/1 , (l ;...; !/ „, i „) = P„_ 1 ( !/1 , il ;...; ! ,„_ 1 , ( „_ 1 ); 

4. j& Ul P l (y 1 M) = \. 

Inasmuch as the distributions P n enable one to compute all the averages of 
the stochastic process [Eq. (j3.4j) ]. they constitute a complete specification of 
it. Conversely, according to a theorem due to Kolmogorov, it is possible to 
prove that the inverse is also true, i.e., that any set of functions obeying the 
above four consistency conditions determines a stochastic process Y(t). 

3.2 Gaussian processes 

A stochastic process is called a Gaussian process, if all its P n are multivariate 
Gaussian distributions (Sec. 12. 3j) . In that case, all cumulants beyond m = 2 
vanish and, recalling that ({Y (h)Y (t 2 ))) = (F(ti)F(t 2 )) - (^(*i)> (Y(t 2 )), 
one sees that a Gaussian process is fully specified by its average (Y(t)) and 
its second moment (Y (ti)Y (t 2 )) ■ Gaussian stochastic processes are often 
used as an approximate description for physical processes, which amounts to 
assuming that the higher-order cumulants are negligible. 
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3.3 Conditional probabilities 

The notion of conditional probability for multivariate distributions can be 
applied to stochastic processes, via the hierarchy of probability distributions 
introduced above. For instance, the conditional probability Pi\i(y2,t2\yi,ti) 
represents the probability that Y takes the value y 2 at t 2 , given that its value 
at t\ "was" y\. It can be constructed as follows: from all sample functions 
Y x (t) of the ensemble representing the stochastic process, select those passing 
through the point y\ at the time t\\ the fraction of this sub- ensemble that goes 
through the gate (7/2, y 2 + ^2) at the time t 2 is precisely -Pi|i (2/2, ^Iz/i, ^i)dy2- 
More generally, one may fix the values of Y at n different times ti, . . . ,t n , 
and ask for the joint probability at m other times t n+ i, . . . , t n+m . This leads 
to the general definition of the conditional probability P m \ n by Bayes' rule: 

p / , , I j. . . j. \ Pn+rniyii ^li • • • i 2/n+ra; tn+m) 

■Lm\n\yn+li ^n+lj • • • i 2/n+m; ^n+m 2/1 ) Z U ■ ■ ■ j yni tn) — 777 T 

P n {yi,ti] ■ . .;y n ,t n ) 
(3-5) 

Note that the right-hand side of this equation is well defined in terms of the 
probability distributions of the hierarchy P n previously introduced. Besides, 
from their consistency conditions it follows the normalisation of the P m \ n . 

3.4 Markov processes 

Among the many possible classes of stochastic processes, there is one that 
merits a special treatment — the so-called Markov processes. 

Recall that, for a stochastic process Y(t), the conditional probability 
Pi\i(y2,t2\yi,ti), is the probability that Y(t 2 ) takes the value y 2 , provided 
Y(ti) has taken the value y±. In terms of this quantity one can express P 2 as 

^Ms/i, h; y 2 , t 2 ) = Piiyi, ti)Pi\i(y 2 , t 2 \yi, h) . (3.6) 

However, to construct the higher-order P n one needs transition 
probabilities P n \ m of higher order, e.g., P3(yi,ti,y 2 ,t 2 ;y 3 ,t 3 ) = 
P 2 (yi, ti, y 2 , t 2 )Pi\ 2 (ys, t 3 \yi, t±] y 2 , t 2 ). A stochastic process is called a Markov 
process, if for any set of n successive times t± < t 2 < • • • < t n , one has 

Pi\n-i(yn,t n \yi,ti;...;y n -i,t n -i) = Pi\i(yn,t n \yn-i,t n -i) ■ (3-7) 

In words: the conditional probability distribution of y n at t n , given the value 
y n -\ at t n -i, is uniquely determined, and is not affected by any knowledge 
of the values at earlier times. 
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A Markov process is therefore fully determined by the two dis- 
tributions Pi(y,t) and P\\\(y' , t'\y, t), from which the entire hierarchy 
PniVii tx'i ■ ■ ■ ] Vn, t n ) can be constructed. For instance, consider ti < t 2 < t 3 ; 
P 3 can be written as 

Psiyi, h; V2, h; y 3 , h) Eq '=' 5 P 2 (yi, h; y 2 , t 2 )P x \ 2 {y 3 , t 3 \y h h; y 2 , t 2 ) 

Eq '=' Pi(Vl, tl\ 2/2, *2)-Pl|l(y3, h\V2, t 2 ) 

Eq.JS) p i (y 1) f 1 )p i|1 (y 2j i 2 |y 1) t 1 )p i|1 (j /gj i 3 |y 2j f 2 ) . (3.8) 

From now on, we shall only deal with Markov processes. Then, the only 
independent conditional probability is Pi\i(y', t'\y, t), so we shall omit the 
subscript 1|1 henceforth and call P\\i(y' ,t'\y,t) the transition probability. 



3.5 Chapman— Kolmogorov equation 

Let us now derive an important identity that must be obeyed by the transition 
probability of any Markov process. On integrating Eq. (13. 8p over y 2 , one 
obtains (ti < t 2 < t 3 ) 

P2(yi,ti,y 3 ,t 3 ) =P 1 (y 1 ,t 1 ) Jdy 2 P(y 2 ,t 2 \y 1 ,t 1 )P(y 3 ,t 3 \y 2 ,t 2 ) , 

where the consistency condition [3] of the hierarchy of distribution functions 
P n has been used to write the left-hand side. Now, on dividing both sides by 
Pi(yi,ti) an d using the special case (I3.6P of Bayes' rule, one gets 



P(ys,ta\yi,h 



<% 2 P(V3, *3 1 2/2 , t 2 )P(y 2 , t 2 \yi, t 



1) , 



(3.9) 



which is called the Chapman-Kolmogorov equation. The time ordering is 
essential: t 2 must lie between t\ and £3 for Eq. (13. 9p to hold. This is required 
in order to the starting Eq. (13.81) being valid, specifically, in order to the sec- 
ond equality there being derivable from the first one by dint of the definition 
( 13.71) of a Markov process. 

Note finally that, on using Eq. (13. 6p one can rewrite the particular case 
Pi (2/2^2) = / dyi P 2 (y 2 , t 2 ; yi, ti) of the relation among the distributions of 
the hierarchy as 



(3.10) 
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This is an additional relation involving the two probability distributions char- 
acterising a Markov process. Reciprocally, any non-negative functions obey- 
ing Eqs. (I3.9P and (13.101) . define a Markov process uniquely. 



3.6 Examples of Markov processes 

Wiener— Levy process. This stochastic process was originally introduced 
in order to describe the behaviour of the position of a free Brownian particle 
in one dimension. On the other hand, it plays a central role in the rigorous 
foundation of the stochastic differential equations. The Wiener-Levy process 
is defined in the range — oo < y < oo and t > through [cf. Eq. (11.51) ] 



P(y2,t 2 \yi,ti) 



_y\_ 

^fhrtl 6XP V 2tx 



: exp 



2vr(t 2 - t h 



(S/2 ~ Vi) 
2(t 2 



ti 



(3.11a) 
{t x <t 2 ) .(3.11b) 



This is a non-stationary (Pi depends on t), Gaussian process. The second- 
order moment is 

(y(ti)y(f 2 )) = min(ti,t 2 ), (3.12) 
Proof: Let us assume t\ <t 2 . Then, from Eq. (13. 4p we have 

(YttJYfo)) = Jd Vl Jdy 2 y 1 y 2 P 2 {y 1 ,t 1 ;y 2 ,t 2 ) 

= JdyiyiPiiyifh) Jdy 2 y 2 P(y 2 ,t 2 \yi,t 1 ) = ydj/iy?Pi(yi,*i) , 



yi by Eq. l !3.11bH h by Eq. fS.llal 

where we have used that t% is the time-dependent variance of Pi. Q.E.D. 

Ornstein— Uhlenbeck process. This stochastic process was constructed 
to describe the behaviour of the velocity of a free Brownian particle in one 
dimension (see Sec. 14.51) . It also describes the position of an overdamped 
particle in an harmonic potential. It is defined by (At = t 2 — 1\ > 0) 

P 1 (y 1 ) = ^=e^{-\yl) , (3.13a) 



P(y2,t 2 \yi,t 1 



2,71 



-2 At 



exp 



(y 2 ~ Vie 
2(1 -e 



-At\2 



-2A^^ 



. (3.13b) 
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The Ornstein-Uhlenbeck process is stationary, Gaussian, and Markovian. 
According to a theorem due to Doob, it is essentially the only process with 
these three properties. Concerning the Gaussian property, it is clear for 
Pi{yx). For P 2 (j/2,*2;2/i,ti) = -Pi (2/1 ) -P (2/2,^2 1 £i) [Eq. flUD], we have 



1 



exp 



y?-2j/iy 2 e A * + 2/2 



2 1 



-2 An 



(3.14) 



This expression can be identified with the bivariate Gaussian distribution 
( I2.23P and the following parameters 



A 



-At 



G' 



1-e 



-2At 



with the particularity that a 2 = 1— A 2 in this case. Therefore, we immediately 
see that the Ornstein-Uhlenbeck process has an exponential auto-correlation 
function (Ytt^Yfo)) = e~ At , since a 2 A/(l - A 2 ) = A in this case@ 

The evolution with time of the distribution P 2 (j/ 2 , * 2 ; £/i, £i), seen as the 
velocity of a Brownian particle, has a clear meaning. At short times the 
velocity is strongly correlated with itself: then A ~ 1 and the distribution 
would be like in the lower right panel of Fig. [5] [with a shrinked variance 
a 2 = (1 — A 2 ) —►()]. As time elapses A decreases and we pass form one panel 
to the previous and, at long times, A ~ and the velocity has lost all memory 
of its value at the initial time due to the collisions and hence P 2 (?/ 2 , Vi, ti) 
is completely uncorrelated. 

Exercise: check by direct integration that the transition probability 
( I3.13bl) obeys the Chapman— Kolgomorov equation (13.91) . 



15 This result can also be obtained by using Eqs. (I3.13|) directly: 



K(h,t 2 ) = (Y(ti)y(ta)) - <y(ti)> (Y(t 2 )) 

= J dyidy 2 yiy2 Pzjyi^ti) yi, fa) 

Pi(yi)P(v2,t 2 \viM) 
3 1 f°° 

dyiyi^=exp(-iyi) / dy 2 y 2 



2tt 



^/2tt(1 - e- 2At ) 



exp 



At\2- 



Oft - yie-n 

2(l-e~ 2At ) 



= e 



-A/ 



dy iy ? exp(-iy2) Q.E.D. 

a V 27T 
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4 The master equation: Kramers— Moyal 
expansion and Fokker— Planck equation 



The Chapman-Kolmogorov equation (13.91) for Markov processes is not of 
much assistance when one searches for solutions of a given problem, because 
it is essentially a property of the solution. However, it can be cast into a 
more useful form — the master equation. 



4.1 The master equation 

The master equation is a differential equation for the transition probabil- 
ity. Accordingly, in order to derive it, one needs first to ascertain how the 
transition probability behaves for short time differences. 

Firstly, on inspecting the Chapman-Kolmogorov equation (13. 9p for equal 
time arguments one finds the natural result 

P{V3, t 3 \y h t) = J dy 2 P(y 3 , t 3 \y 2 , t)P(y 2 , t\y x , t) P(y 2 , t\y lt t) = 8{y 2 -y x ) 

which is the zeroth-order term in the short-time behaviour of P(y',t'\y,t). 
Keeping this in mind one adopts the following expression for the short-time 
transition probability: 



P(y 2 , t + At\ Vl , t) = 5(y 2 - Vl ) [1 - a<°> (y h t) At] + W t (y 2 \ yi )At + 0[(At) 



(4-1) 

where Wt(y 2 \yi) is interpreted as the transition probability per unit time from 
yi to y 2 at time t. Then, the coefficient 1 — a^(y%, t)At is to be interpreted 
as the probability that no "transition" takes place during At. Indeed, from 
the normalisation of P (y 2 , t 2 \yi, ti) one has: 

1 = Jdy 2 P{y 2 ,t+ At\ yi ,t) ~ l-a (0) (2/i,t)At + Jdy 2 W t {y 2 \y x ) At . 

Therefore, to first order in At, one get 3^1 

a^(y u t) = Jdy 2 W t {y 2 \ yi ) , (4.2) 

16 Thc reason for the notation will become clear below. 
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which substantiates the interpretation mentioned: a^°'(yi,t)At is the total 
probability of escape from yi in the time interval (t, t + At) and, thus, 1 — 
a^(yi, t)At is the probability that no transition takes place during this time. 

Now we can derive the differential equation for the transition probability 
from the Chapman-Kolmogorov equation (13.91) . Insertion of the above short- 
time expression for the transition probability in into it yields 

%3-W2)[l-a C0) (»2,t2)At]+W t2 (i/ 3 |ite)Af 

A 



P{y 3 , t 2 + At\ Vl , U) = J dy 2 P(y 3 , t 2 + At\y 2 , t 2 ) P{y 2 , t 2 \y u h) 

~ [l-a^{y^t 2 )^t]P{yzMvi,h) 

+ At Jdy 2 W t2 (y 3 \y 2 )P{y 2 ,t 2 \y 1 ,t 1 ) . 

Next, on using Eq. (14.21) to write a^°'(y 3 ,t 2 ) in terms of W t2 (y 2 \y3), one has 
^ [P(y 3 , h + At\ Vl , h) - P(y 3 , t 2 \y h h)} 

dy 2 [W t2 {y 3 \y 2 )P{y2,t 2 \yi,t 1 ) - W t2 {y 2 \y 3 )P{y 3 ,t 2 \yiM)\ 



which in the limit At — ► yields, after some changes in notation (yi,t\ 
yo,to, 2/2, ^2 — ► y',t, and y 3 — ► y), the master equation 




(4.3) 

which is an integro-differential equation. 

The master equation is a differential form of the Chapman-Kolmogorov 
equation (and sometimes it is referred to as such). Therefore, it is an equa- 
tion for the transition probability P(y, t\yo, to), but not for Px(y, t). However, 
an equation for Pi(y,t) can be obtained by using the concept of "extraction 
of a sub-ensemble". Suppose that Y(t) is a stationary Markov process char- 
acterised by Pi(y) and P(y,t\y ,t ). Let us define a new, non-stationary 
Markov process Y*(t) for t > t by setting 

P^yuh) = F(yi,ti|2/ ,to) , (4.4a) 

P*(v2,h\yi,ti) = PfaMvuti) ■ ( 4 -4b) 
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This is a sub-ensemble of Y(t) characterised by taking the sharp value y 
at t , since P^(yi,t ) = S(yi — y ). More generally, one may extract a sub- 
ensemble in which at a given time to the values of Y*(to) are distributed 
according to a given probability distribution p(yo)'- 

Pi{yi,h) = fdyo P(yi,h |yo, t )p(y ) , (4.5) 



and P*(y2,t 2 \yi,ti) as in Eq. (14.4bj) . Physically, the extraction of a sub- 



ensemble means that one "prepares" the system in a certain non-equilibrium 
state at to. 

By construction, the above P*(yi,ti) obey the same differential equa- 
tion as the transition probability (with respect its first pair of arguments), 
that is, P*(yi,ti) obeys the master equation. Consequently, we may write, 
suppressing unessential indices, 



dy' [W(y\y')P(y', t) - W(y'\y)P(y, t)) . (4.6) 



dP(y,t) 
dt 

If the range of Y is a discrete set of states labelled with n, the equation 
reduces to 

= Wnn<Pn>{t) ~ W n , nPn (t)} . (4.7) 
n' 

In this form the meaning becomes clear: the master equation is a balance 
(gain-loss) equation for the probability of each state. The first term is the 
"gain" due to "transitions" from other "states" n' to n, and the second term 
is the "loss" due to "transitions" into other configurations. Remember that 
Wn'n > and that the term with n = n' does not contribute. 

Owing to W(y\y')At is the transition probability in a short time interval 
At, it can be computed, for the system under study, by means of any available 
method valid for short times, e.g., by Dirac's time-dependent perturbation 
theory leading to the "golden rule". Then, the master equation serves to 
determine the time evolution of the system over long time periods, at the 
expense of assuming the Markov property. 

The master equation can readily be extended to the case of a multi- 
component Markov process Yi(t), i = 1, 2, . . ., N, on noting that the Chapman- 
Kolmogorov equation (13. 9p is valid as it stands by merely replacing y by 
y = (z/ij ' ' " ; 2/jv) • Then, manipulations similar as those leading to Eq. (14.61) 
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yield the multivariate counterpart of the master equation 




(4.8) 



Example: the decay process. Let us consider an typical example of 
master equation describing a decay process, in which p n (t) determines the 
probability of having at time t, n surviving "emitters" (radioactive nuclei, 
excited atoms emitting photons, etc.). The transition probability in a short 
interval is 

{0 for n> n' 
772' At for n = n' -1 
O(At) 2 for n< n' - 1 

That is, there are not transitions to a state with more emitters (they can 
only decay; reabsortion is negligible), and the decay probability of more that 
one decay in At is of higher order in At. The decay parameter 7 can be 
computed with quantum mechanical techniques. The corresponding master 
equation is Eq. (14.71) with W n%n i = jn'S n n '^i 



and hence 



dt 

dt 



,n+l Pn+1 (t) ~ Wn-i, n p n (t) ■ 



7(n + 1) p n +i(t) - jnp n (t) 



(4.9) 



Without finding the complete solution for p n (t), we can derive the equa- 
tion for the average number of surviving emitters (N) (t) = J2^=o n Pn(t) 



k=n+l 



^n(dp n /dt) = 7^n(ra + l)p n+ i -7^n 2 p n 

n=0 



n=0 



n=0 



7 



<T,[(V-i)k- f\ 



Pk 




k=l 



k=l 



Therefore the differential equation for (N) and its solution are: 



_ ( N ) = -7 (N) 



(N) (t) = n Q e 



-it 



(4.10) 
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4.2 The Kramers— Moyal expansion and the Fokker— 
Planck equation 

The Kramers-Moyal expansion of the master equation casts this integro-dif- 
ferential equation into the form of a differential equation of infinite order. It 
is therefore not easier to handle but, under certain conditions, one may break 
off after a suitable number of terms. When this is done after the second-order 
terms one gets a partial differential equation of second order for P(y, t) called 
the Fokker-Planck equation. 

Let us first express the transition probability W as a function of the size 
r of the jump from one configuration y' to another one y, and of the starting 
point y'\ 

W(y\y')=W(y';r) , r = y - y' . (4.11) 
The master equation (14.61) then reads, 

£ = JdrW(y-r;r)P(y-r,t)-P(y,t)JdrW(y;-r), (4.12) 

where the sign change associated with the change of variables y' — > r = 
y — y' i is absorbed in the boundaries (integration limits), by considering a 
symmetrical integration interval extending from — oo to oo: 



dy'f(y') 



py—oo /' — oc poo 

- drf(y-r) = - drf(y-r)= drf{y 

J J/+OO J oo J —oo 



Moreover, since finite integration limits would incorporate an additional de- 
pendence on y, we shall restrict our attention to problems to which the 
boundary is irrelevant. 

Let us now assume that the changes on y occur via small jumps, i.e., that 
W(y'\ r) is a sharply peaked function of r but varies slowly enough with y' . 
A second assumption is that P(y, t) itself also varies slowly with y. It is then 
possible to deal with the shift from y to y — r in the first integral in Eq. 
(I4.12p by means of a Taylor expansion: 



dP(y,t) 

at 



drW(y;r)P(y,t) + J2 



m=l 



drr r 



d" 



dy r 



■[W(y,r)P(y,t)\ 



-P(y,t) JdrW{y;-r) 



E 

m=l 



-1 



gr, 



m am 



ml dy r 



drr m W(y; r) 



P(y,t) 
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where we have used that the first and third terms on the first right-hand 
side cancel each other0 Note that the dependence of W(y; r) on its second 
argument r is fully kept; an expansion with respect to it, is not useful as W 
varies rapidly with r. Finally, on introducing the jump moments 



,( m )i 



(y,t) = drr m W(y;r) 



(4.13) 



one gets the Kramer s-Moyal expansion of the master equation: 



dP(y,t) 
dt 



'-l) m d n 



4 m\ dy r - 



[a^( y ,t)P( y ,t)] . 



(4.14) 



Formally, Eq. (14.141) is identical with the master equation and is therefore 
not easier to deal with, but it suggest that one may break off after a suitable 
number of terms. For instance, there could be situations where, for m > 2, 
a( m \y,t) is identically zero or negligible. In this case one is left with 



dP(y,t) 
dt 



(2) 



(y,t)P(y,t)] 



(4.15) 



which is the celebrated Fokker-Planck equation. The first term is called the 
drift or transport term and the second one the diffusion term, while a^(y, t) 
and a^ 2 '(y,t) are the drift and diffusion "coefficients". 

It is worth recalling that, being derived from the master equation, the 
Kramers-Moyal expansion, and the Fokker-Planck equation as a special case 
of it, involve the transition probability P(y, t\yo, to) of the Markov stochastic 
process, not its one-time probability distribution P\(y, t). However, they also 
apply to the P*{y, t) of every subprocess that can be extracted from a Markov 
stochastic process by imposing an initial condition [see Eqs. (14.41) and (14.51) 



4.3 The jump moments 

The transition probability per unit time W(y'\y) enters in the definition 
( I4.13P of the jump moments. Therefore, in order to calculate a^ m \y,t), we 

17 This can be shown upon interchanging — r by r and absorbing the sign change in the 
integration limits, as discussed above. 
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must use the relation (14.1 j) between VK(y'|y) and the transition probability 
for short time differences. 

Firstly, from Eq. (14.111) one sees that W(y'; r) = W(y\y') with y = y' + r. 
Accordingly, one can write 



W(y;r) = W(y'\y) , 



y =y + r 



On inserting this expression in Eq. (14.131) one can write the jump moments 
JTsI 



(4.16) 




In order to calculate the jumps moments we introduce the quantity 
A^\ y] T,t)= [dy'(y'-y) m P(y',t + r\y,t), (m > 1) , 



which is the average of [Y(t + r) — Y(t)] m with sharp initial value Y(t) = y 
(conditional average). Then, by using the short-time transition probability 
(14.11) . one can write 



A {m \y;r,t) = jdy'(y'-yr{5(y'-y)[l-a^(y,t)r] + W(y'\y)T + 0(r 2 )} 
= Tjdy'(y'-yrW(y'\y) + 0(r 2 ) 



a^(y,t)r + 0(r 2 ) 



(m > 1) 



where the integral involving the first term in the short-time transition prob- 
ability vanishes due to the presence of the Dirac delta. Therefore, one can 
calculate the jump moments from the derivatives of the conditional averages 
as follows 

d 



(m) 



(y,t)= — A^(y,T,t) 



r=0 



Finally, on writing 

A ( ~ m \y;At,t) = Jdy' (y'-y) m P(y' ,t+At\y,t) = ([Y(t + At) - Y(t)}' 



Y{t)=y 



18 This equation makes clear the notation employed. The quantity a^°\y,t) [Eq. (14. 2| 
which was introduced in Eq. (|4.1|) . is indeed the m = jump moment. 
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one can alternatively express the jump moments as 



iT„Ai< [y(t + At) - y(t)l ' 



Y(t)=y 



(4.17) 



In Sec. below, which is devoted to the Langevin equation, we shall calcu- 
late the corresponding jump moments in terms of the short-time conditional 
averages by means of this formula. 

4.4 Expressions for the multivariate case 

The above formulae can be extended to the case of a multi-component 
Markov process Yi(t), i = 1, 2, . . ., N. Concerning the Kramers-Moyal ex- 
pansion one only needs to use the multivariate Taylor expansion to get 



dP 
~dt 



y^-y 



(m) 



(4.18) 



while the Fokker-Planck equation is then given by 



dP 



< >r s. — ■> d ( i ) , j \ - 



dyi 



Q2 



2 ^ d yi d yj 



a%\y,t)P 



(4.19) 



In these equations, the jump moments are given by the natural generalisation 
of Eq. (14.161) , namely 

41iJy> *) = / & (y'n -v*)'~ (vL - yjjw(y'\y) , (4.20) 

and can be calculated by means of the corresponding generalisation of Eq. 



1 / 

(y > = in xt \ II ^ (* + At ) - ^ 



Y k (t)=y k 



(4.21) 



that is, by means of the derivative of the corresponding conditional average. 
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4.5 Examples of Fokker— Planck equations 

Diffusion equation for the position. In Einstein's explanation of Brow- 
nian motion he arrived at an equation of the form [see Eq. (jl.4p 



(4.22) 




Comparing with the Fokker-Planck equation (I4.15p . we see that in this case 
a^'(x, t) = 0, since no forces act on the particle and hence the net drift 
is zero. Similarly a^(x,t) = 2D, which is independent of space and time. 
This is because the properties of the surrounding medium are homogeneous 
[otherwise D = D(x)]. The solution of this equation for P(x,t = 0) = S(x) 
was Eq. (II .5p . which corresponds to the Wiener-Levy process (13.111) . 

This equation is a special case of the Smoluchowski equation for a parti- 
cle with large damping coefficient 7 (overdamped particle), the special case 
corresponding to no forces acting on the particle. 




Diffusion equation in phase space (x,v). The true diffusion equation 
of a free Brownian particle is 



(4.23) 



This equation is the no potential limit of the Klein-Kramers equation for 
a particle with an arbitrary damping coefficient 7. From this equation one 
can obtain the diffusion equation (14.221) using singular perturbation theory, 
as the leading term in a expansion in powers of I/7. Alternatively, we shall 
give a proof of this in the context of the Langevin equations corresponding 
to these Fokker-Planck equations!^! 

We have stated without proof that the Ornstein-Uhlenbeck process de- 
scribes the time evolution of the transition probability of the velocity of 
a free Brownian particle. We shall demonstrate this, by solving the equa- 
tion for the marginal distribution for v obtained from (I4.23p . The marginal 



19 We shall see that the Langevin equation m'x = —mjx + £(t) [Eq. (|1 -8|t ] leads to Eq. 
(|4.23p , while the overdamped approximation mj x ~ £(i) corresponds to Eq. (14.221) . 
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probability is Py(v,t) = fdxP(x, v, t). Integrating Eq. (14.231) over x, using 
/ dx d x P(x, v, i) = 0, since P(x = ±00, v, t) = 0, we find 

dP v ( d k B T d 2 \ n , t A . 

7 [^-v + -?-—)P v . (4.24) 



dt \dv m dv 2 r 

We will see that this equation also describes the position of an overdamped 
particle in an harmonic potential. Thus, let us find the solution of the generic 
equation 

' rd t P = d y {yP) + Dd 2 P . I (4.25) 



Introducing the characteristic function (12. 2p (so we are solving by the Fourier 
transform method) 



G(k,t) = J dye ik yp(y,t) , P(y,t) = Lj dke~ ik VG(k, 



t) 



the second order partial differential equation (14.251) transforms into a first 
order one 

rd t G + kd k G = -Dk 2 G , (4.26) 

which can be solved by means of the method of characteristics!^! 
In this case the subsidiary system is 

dt dA; dG 



r k Dk 2 G ' 
Two integrals are easily obtained considering the systems t, k and k, G: 

— = ^ — > k = a e l l T — > u = ke~^ T = a 
t k 

-Dkdk = dG/G -> -\Dk 2 = lnG + c -> v = e'z^G = b 



20 In brief, if we have a differential equation of the form 

P d J f + Q d 7 f = R, 
ox ay 

and u(x, y, f) = a and v(x, y, /) = b are two solutions of the subsidiary system 

da; dy df 
~P = ~Q = ~R ' 

the general solution of the original equation is an arbitrary function of u and v, h(u, v) = 0. 
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Then, the solution h(u, v) = can be solved for v as v = <fi{u) with still an 
arbitrary function 0, leading the desired general solution of Eq. (14.261) 



G 



-■kDk 2 
e 2 



>(fce-^) , (4.27) 

by means of the methods of characteristics. 

The solution for sharp initial values P(y,t = 0) = S(y — y ) leads to 
G(k, t = 0) = exp(iky ), from which we get the functional form of 0: (j>(k) = 
exp(iky + ^Dk 2 ). Therefore, one finally obtains for G(k,t) 

G{k, t) = exp [iy e- t/T k - \D{1 - e- 2t ' T ) k 2 } , (4.28) 

which is the characteristic function of a Gaussian distribution [see Eq. (12.91) ]. 
with Hi = y§e~ t l T and o 2 = D(l — e _2 '/ r ). Therefore, the probability distri- 
bution solving Eq. (14.251) is 



P{y,t\y ,o) 



'1tiD(\ - e- 2t / T ) 



exp 



(y-ype t/T ) 2 

2D(1 -e~ 2t l T ) 



(4.29) 



which, as stated, is the transition probability of the Ornstein-Uhlenbeck 
process [Eq. ( l3~T3bl) ]. Q.E.D. 

Note finally that the parameters for the original equation for Py [Eq. 
Pi) ], are simply \i x = v e~ t/T and a 2 = (k B T/m)(l - e" 2 */ T ). Thus, at 
long times we have Py oc exp(— \mv 2 jk^T) which is simply the statistical 
mechanical equilibrium Boltzmann distribution for free particles. 



Diffusion equation for a dipole. The diffusion equation for a dipole 
moment p in an electric field E is (neglecting inertial effects) 



dP 



1 d 



dt sin $ di!) 



OP 

k^T— — h pE sin-^P 
ov 



(4.30) 



This equation was introduced by Debye in the 1920's, and constitutes the 
first example of rotational Brownian motion. ( is the viscosity coefficient 
(the equivalent to 7 in translational problems). It is easily seen that P oc 
exp(pE cos$/ k B T) is the stationary solution of Eq. (14.301) . which leads to 
the famous result for the average dipole moment of an assembly of dipoles, 
(cos??) = cotha — 1/a with a = pEjk^T and to Curie's paramagnetic law 
at low fields. However, Eq. ( I4.30p also governs non-equilibrium situations, 
and in particular the time evolution between different equilibrium states. 
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5 The Langevin equation 



5.1 Langevin equation for one variable 

The Langevin equation for one variable is a "differential equation" of the 
form [cf. Eq. (Ofr 



^L=A(y,t) + B(y,t)m 



(5.1) 



where is a given stochastic process. The choice for £(t) that renders 

a Markov process is that of the Langevin "process" (white noise), which is 

Gaussian and its statistical properties are 



(mm) 



o, 

2L>5(t 1 - t 2 



(5.2a) 
(5.2b) 



Since Eq. ( 15. ip is a first-order differential equation, for each sample function 
(realisation) of £(t), it determines y(t) uniquely when y(t ) is given. In ad- 
dition, the values of the fluctuating term at different times are statistically 
independent, due to the delta-correlated nature of £(t). Therefore, the val- 
ues of £(t) at previous times, say t' < t , cannot influence the conditional 
probabilities at times t > to- From these arguments it follows the Markovian 
character of the solution of the Langevin equation (15. ip . 

The terms A(y, t) and B(y, t)^(t) are often referred to as the drift (trans- 
port) and diffusion terms, respectively. Due to the presence of £ (t) , Eq. (15.11) 
is a stochastic differential equation, that is, a differential equation comprising 
random terms with given stochastic properties. To solve a Langevin equation 
then means to determine the statistical properties of the process y(t). 

Finally, the higher-order moments of £(t) are obtained from the second 
order ones (15.21) . by assuming relations like those of the multivariate Gaussian 
case, i.e., all odd moments of £(t) vanish and, e.g., 

(mmmm) = m 2 [ (mm) (mm) 

+ (mm) (mm) 
+ (mm)(mm)] ■ (5.3) 



21 Hereafter, we use the same symbol for the stochastic process Y(t) and its realisations 

y(t). 
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To check this result, we shall demonstrate a slightly more general result 
known as Novikov theorem. 



5.1.1 Novikov theorem and Wick formula 



Novikov theorem states that for a multivariate Gaussian distribution (Sec. 
12.31) with zero mean 



P(x) 



detA 
(2ir) n 



exp (— |x • A ■ x 



the averages of the type (xj/(x)), can be obtained as 




(5.4) 



(5.5) 



Applying this result to /(x) = XjXkX? and using dxi/dx m = 5i m , we have 

t\X iX j X foX ^ {•Ei'^m 

df/dx m 

Therefore, using the Kronecker's delta to do the sum, we get Wick's formula 

(xiXjXkXt) = (xiXj) (x k x e ) + (xiX k ) (xjXe) + {xiXi) (xjX k ) . (5.6) 

Equation ( 15.31) then follows because £(£) is assumed to be a Gaussian pro- 
cess, which by definition means that the n times probability distribution 
Pa (£i, t\\ . . . ; £ n , t n ) is a multivariate Gaussian distribution. 

Proof of Novikov theorem. We shall demonstrate this theorem in three 
simple steps: 

(1) If we denote by E{x) = | . x^A^Xj minus the exponent in the Gaussian 
distribution (15.41) . we have 

^ XiAijXj 2 ^ ^ ^imAjjXj -\- Xj^4.jj5j m ) 



dx r 



2 dx r 



[A^ sym.] 



2 A m jXj + 2 XjA 
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(2) Using the definition of average for (xj/(x)), inserting the above expression 
for Xi and integrating by parts, we have 



(xj(x)> = C jdxxJ(x)e- E 

= C^^" 1 )™ /dx/(x 

m J 





(3) Finally, we demonstrate that (XiXj) = (A (a particular case of 
the result {(XiXj)) = given without proof in Sec. 12.31) . Indeed, 

using the above result for / = Xj and dxj/dx m = 8j m , we have (xiXj) = 
^ m (y4 _1 ) im 5j m = (A~ l )ij. Insertion of this in the above result completes the 
proof of Novikov's theorem. 

5.2 The Kramers— Moyal coefficients for the Langevin 
equation 

Since the solution of the Langevin equation is a Markov process, it obeys a 
master equation, which may be written in the Kramers-Moyal form (14. 14ft . 
Let us calculate the successive coefficients (14.171) occurring in that expan- 
sion. We first cast the differential equation (15. ip into the form of an integral 
equation 



/t+At pt + lAt 

dt 1 A[y(t 1 ),t 1 }+ J dti£[y(ti),ti]f(ti) , (5.7) 

where y stands for the initial value y(t). On expanding according to 

>%(ti),ti] = A{y,t 1 ) + A'{y,t 1 )[y{t 1 )~y] + --- , 
B[y{t x ),t x ] = B(y,t 1 )+B'(y,t 1 )[y(t 1 )-y) + --- , 

where the prime denotes partial derivative with respect to y evaluated at the 
initial point: 



ft+At 



At, \ OA 
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one gets 



-t+At 

!i(t-\ St)-/, = I dt x A{yM) 

-t+At 



/t+i\t 
dt x A'(y,t 1 )[y(t 1 )-y) + 



/t+At 
dhBiyMWi) 

rt+At 

+ / dt 1 B'(y,t 1 )[y(t 1 )-y}£(t 1 ) + --- . (5.8) 



For y(ti) — y in the above integrands we iterate Eq. (15. 8p to get 

/t+At 
dt x A{yM) 

t+At i-ti 



+ 
+ 



5 /-tl 

dtiA'fati) / dt 2 A(y,t 2 ) 

t rti 

dt x A{yM) J dt 2 B(y,t 2 W 2 ) + 



rt+At 

+ / dtiB(y,ti)e(*i) 



dttB'iy^ih) J dt 2 A(y,t 2 ) 

/t+At rt! 
dh B'(y, h^ih) J dt 2 B(y, t 2 )£(t 2 ) + • • ■ (5.9) 

If we take the average of this equation for fixed y = y(t), by using the 
statistical properties ( \5.2L we obtain the conditional average required to get 

/t+At rt+At rti 

dt l A(y,t 1 ) + J te x A\yM)] dt 2 A(y,t 2 ) 

/t+At pti 
dtiB'{y,ti)J dt 2 B{y,t 2 )5{t 2 -t 1 ) + -- - . 
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Next, on using for the Dirac delta the result / dt8(t — t )f(t) = |/(to)j we 

Jt 

obtain 

f ti 

J dt 2 B{y, t 2 )5(t 2 - t x ) = \B{y, h) . (5.10) 

Finally, on considering that = ^im ^ (y(t + At) — y) \ y (t)= y , for the cal- 
culation of which only terms through order At need to be retained, one finally 
gets 

a^(y,t)=A(y,t) + DB(y,t)^^-. 

Other integrals not written down in the above formulae do not contribute 
in the limit At — > 0. This can be seen as follows: each Langevin fluctuating 
term on the right-hand side of Eq. (15.91) . is accompanied by an integral. The 
lowest-order terms are written in that expression, whereas higher-order terms 
can be of two types: (i) Integrals of the form of, e.g., 

( / dh ■ / dt 2 ■ ..£(t 2 ) / dt 3 ■ ■ .£(* 3 ) / dt 4 ■ ■ .£(U) ) , 

\Jt Jt Jt Jt I 

which can only give a contribution proportional to (At) 2 , as it is seen by 
using the splitting of (£ (ii) £ (£2) £(^3) m sum of products of the form 
(C(tk)^(ti)) [Eq. f)5.3p ]. (ii) Integrals containing no Langevin 
terms, which are proportional to (At) n , where n is the number of simple 
integrals. Both types of terms clearly vanish when dividing by At and tak- 
ing the limit At — > 0. 

On using the same type of arguments to identify some vanishing integrals 
one can compute the second coefficient in the Kramers-Moyal expansion, 
a( 2 ) = jim i ([y(t + At) - yf) \ y(t)=y , obtaining 

^ rt+At rt+At 

a (2) (y,t) = }im o — J dtxBfatx) J dt 2 B(y, t 2 )2D(5(t 1 -t 2 ) = 2DB 2 (y,t) 
whereas all the coefficients vanish for m > 3. Thus, on collecting all 
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these results one can finally write 



a 



A(y,t) + DB(y,t) 



dB(y,t) 
dy 



(5.11) 



a 



a {2 \y,t) 
i {m) (y,t) 



2DB 2 (y,t), 







for m > 3 . 



5.3 Fokker— Planck equation for the Langevin equation 

From Eq. (15.111) it follows that, for the Markov stochastic process deter- 
mined by the Langevin equation (15. ip with Gaussian ^-correlated £(t), the 
Kramers-Moyal expansion includes up to second-order terms. Therefore, 
the distribution of probability obeys a Fokker-Planck equation [Eq. (14.151) ]. 
which in terms of the above jump moments is explicitly given by 



term, DB(y,t)B'(y,t), which is called the noise-induced drift. This equation 
is very important, since it allows one to construct the Fokker-Planck equation 
directly in terms of the coefficients appearing in the equation of motion. In 
some cases, it can even be done by simply inspection of that equation. 

5.3.1 Multivariate case 

The stochastic differential (Langevin) equation for a multi-component pro- 
cess y = (yi, • • • , yx) has the form 





(5.13) 



k 
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where the are N L white- noise termso The statistical properties of the 
are 



(&(*)> = 0, 
(Uti)Uh)) = 2D5 u 5{h-t 2 ) . 



(5.14a) 
(5.14b) 



Again, the higher-order moments are obtained from these ones, on assuming 
relations like those of the (multivariate) Gaussian case. 

The successive coefficients (I4.2ip occurring in the Kramers-Moyal expan- 
sion (14.181) can be calculated by using arguments entirely analogous to those 
employed above to identify some vanishing integrals. On doing so, one gets 
the following generalisation of Eqs. (15.111) in the multivariate case: 



,(!), 



y,t) = Ai(y,t)+Dj2 B jk (y,t) 



jk 



dB ik (y,t) 



,(2) 



(y,t) = 2Dj2 B ik(y,t)B jk (y,t) 



(5.15) 



(m) 



o 



for m > 3 . 



Again, for the Markov stochastic process defined by the set (15.131) of Langevin 
equations, the Kramers-Moyal expansion of the master equation includes up 
to second-order terms, so that the probability distribution obeys a Fokker- 
Planck equation 



i v L jk 



p 



9 2 



dyidyj 



B ik (y,t)B jk (y,t) 



k 



P 



(5.16) 



which is entirely determined by the coefficients of the Langevin equation. 



22 The number of Langevin sources, N-^, does not need to be equal to the number of 
equations. For example, the sum in k in Eq. (|5 . 13|) can even have one term, ./Vl = 1 — the 
case of "scalar noise" . 
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5.4 Examples of Langevin equations and derivation of 
their Fokker— Planck equations 

5.4.1 Diffusion in phase-space: Klein— Kramers equation 

Let us consider the following generalisation of the original Langevin equation 
(11.81) . in order to account for the presence of an external potential U(x,t) 
(e.g., gravity in the Brownian motion problem) 



d 2 x dx dU . 
m^-rr = —my — h m£{t) . 



dt 2 



dt dx 



(5.17) 



This is simply Newton equation augmented by the fluctuating force [for con- 
venience we have extracted m from £(£)]. 

Let us divide by the mass, introduce V = U/m and the notation V = 
dV/dx, and write (I5.17P as a pair of first-order differential equations 



dx 

dt 
dv 

dt 



(5.18) 
(5.19) 



Then, comparing with the multivariate Langevin equation (15.13j) . we identify 
£ x (t) = and £ v (t) = as well as 



A, 
A, 



-(yv + V) 



P*xx — P*xv — 
Pyx — Pyv 1 



Inserting these results in the general Fokker-Planck equation (15.161) . one gets 
(note djBik = 0) 



dP 
~dt 



d d d 2 

._„__ H7t , + ni+fl _ 



P . 



Gathering the Hamiltonian terms and identifying D/y — k^T/m, we finally 
find the famous Klein-Kramers equation 



(5.20) 
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The result D/^f = k^T/m comes from inserting the Boltzmann distribu- 
tion P oc exp[— (|mu 2 + mV)/k-QT] and finding the conditions for it to 
be a stationary solution. This is equivalent to P. Langevin recourse to the 
equipartition theorem (Sec. 11.1.21) to find (mv 2 ) = k^T. Note finally that 
in the absence of potential, the Klein-Kramers equation leads to the equa- 
tion for free diffusion (I4.23p . with solution (for the marginal distribution 
P v (v,t) = fdxP(x,v,t)) given by the Ornstein-Uhlenbeck process (14.291) . 



5.4.2 Overdamped particle: Smoluchowski equation 

Let us consider the overdamped limit of the Newton-Langevin equation 
(EUD 



(5.21) 



Comparing with the univariate Langevin equation (15.11) . we identify 

A = -V'h , B = I/7 . 




Inserting these results in the Fokker-Planck equation (I5.12p . one gets (putting 
again D/j — k-^I '/m) the Smoluchowski equation 



(5.22) 



The result D/'j = k^T/m can also be obtained on inserting the marginal 
Boltzmann distribution P$ oc exp[— mV{x)/k-BT] and finding the conditions 
for it to be a stationary solution. 

In the absence of potential, the Smoluchowski equation leads to the equa- 
tion for free diffusion (14.221) or Einstein's Eq. (I1.4p . with solution given by the 
Wiener-Levy process [Eqs. (II. 5p or (13. lip ]. Note also that in an harmonic 
potential V(x) = ^uj 2 x 2 , the equation is equivalent to Eq. (I4.25p . with param- 
eters t = 7/c^o> D = kftT/muj 2 , whose solution is the Ornstein-Uhlenbeck 
process (I4.29p . Therefore we can at once write for the overdamped harmonic 
oscillator 



P(x,t) 



mult. 



2nk B T(l - e 



-2t/r' 



exp 



muj 2 (x — x e t l T ) 2 
' 2fc R T(l - e~ 2t l T ) 



T 



2 

• 



7M 

(5.23) 
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Thus, at long times we have P oc exp(— ^mujQX 2 / ksT) which is simply the 
statistical mechanical equilibrium Boltzmann distribution for the harmonic 
oscillator. In addition Eq. f!5.23|) tells us how the relaxation to the equilibrium 
state proceeds. 



5.4.3 Equations for a classical spin (dipole) 




For classical spins, the Langevin equation is the stochastic Landau-Lifshitz 
equation, which for instance describes the dynamics of the magnetic moment 
of a magnetic nanoparticle. Written in simplified units (fields in frequency 
units), it reads 



(5.24) 



Here, B = —dTi/dsis the effective field associated with the Hamiltonian 
of the spin 7i(s) (the equivalent to F = —dll/dx in mechanical problems), 
and the double vector product is the damping term, which rotates s towards 
the potential minima (preserving its length). The stochastic properties of the 
components of £(t) are the usual ones [Eq. f)5.14p ]. but £(£) is now interpreted 
as a fluctuating field. Finally, the damping coefficient A measures the relative 
importance of the relaxation and precession terms. 

The stochastic Landau-Lifshitz equation ( 15.241) . can be cast into the form 
of the general system of Langevin equations (I5.13p . by identifying 



Ai = tjjk SjB k + A y^(g 2 ^ fc - SjS k )B k 



B{ k ^ ^'ijk^j 



(5.25) 



(5.26) 



where is the antisymmetrical unit tensor of rank three (Levi-Civita sym- 
bol^ and we have expanded the triple vector products — s A (sAB) by using 
the rule a A (6 A c) = 6(2 ■ c) - c(a-b) ("B AC-CAB" rule). 

To calculate the noise-induced drift coefficient of the Fokker-Planck equa- 
tion [the term accompanying Ai in Eq. f)5.16p ] we need the derivative of the 



23 This tensor is defined as the tensor antisymmetrical in all three indices with e xyz = 1. 
Therefore, one can write the vector product of A and B as (^4 A B) i = Yljk e ijkAjBk- In 
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diffusion "coefficient" : 



dB 



ik 



ds.j 



-ijk 



(5.29) 



On multiplying by Eq. (15.261) for summing, and using the second con- 
traction property (15.281) for the e^-fc, one finds 



jk 3 £ jk 




Let us compute now the coefficient in the diffusion term Y^ fc B ik Bj k : 



E = E ( E ( E = E *». E 

k k r s r,s k 



where we have employed the contraction rule (15.271) . 

Therefore the Langevin equation associated to the stochastic Landau- 
Lifshitz equation (I5.24p . reads 

dP d f r 11 

-^T = Y e ijkSjB k + X^^dik - SiS k )B k - 2Dsi Pj 



o 2 



dsidsj 



S S'lj Si Sj ] P 



Taking the Sj-derivative in the last term by using ^2jdj(s 2 5ij — s^Sj) - 
Y^j( 2s jSij-$ij s j- s idjj) = 2si-Si-3si, we obtain D £E d i [-2s i P+Y Jj (s 2 5 ij - 



addition, one has the useful contraction property 



^ ^ ^ijk^i' y k $ii'$jj ! $ij'$ji' 
k 

^ ^ ^-ijk^-i'jk *2Sa' j ^ \ ^ijk tjjk 6 . 

jk ijk 

where the last two are obtained by repeated contraction of the first one. 



(5.27) 
(5.28) 
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SiSj)djP]. The first term cancels —2Dsi in the drift and the 2nd can be com- 
bined with *Yli k (s 2 5ik — SiSk)Bk- Finally, returning to a vector notation 



dP 
~dt 



d_ 

ds 



s A B - Xs A 



s A IB — k n T 



0_ 

ds 



■P 



(5.30) 



where (d/ds) ■ J = J2i(^^i/^ s i) (divergence) and, by analogy with the me- 
chanical problems, we have set D = Xk-gT. This equation can be seen as the 
rotational counterpart of the Klein-Kramers equation. 

The electric case corresponds to the precession term dominated by the 
damping term (a sort of Smoluchowski equation) 



dP 
~dt 



X 







p A 



pA ( E — k B T 







>P . 



(5.31) 



Then, introducing spherical coordinates ip) and assuming the Hamiltonian 
to be axially symmetric, the equation above reduces to 

d 



1 dP 
X~dt 



1 



sin •& dd 



on dp\ 



(5.32) 



This equation corresponds to Eq. ( 14.30P by identifying 1/A — > ( and H = 
—pE cos i9, which is the Hamiltonian of the dipole in an external field. 

Let us solve Eq. (14.301) for a time dependent field E(t) with E(t < 0) = 
AE, while E(t > 0) = 0. That is, in the distant past the system was 
equilibrated in the presence of a small field AE which at t = is removed, 
and we seek for the time evolution of P($, t). For t < it is easily seen 
that Pq = N exp(pAE cos $/k B T) is the stationary solution of Eq. (I4.30p . 
since dPo/dfl = —{pAE said / k-QT)P . Then, introducing the notation a = 
pAE/k B T, and using that q < 1, we have 

P (tf) ~iV(l + a;costf) , (t < 0) . (5.33) 

For t > the field is removed and for the solution we use the ansatz 



P(&, t) = N[l + ag(t) costf] , (t > 0) 



(5.34) 



with g(t) a function to be determined by inserting this P in the Fokker- 
Planck equation (jODjl with E(t > 0) = 0: 



C^a-TT cos ^ = ~k B TNa g— 

at sm v ov 



° (sin 2 ^) 
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Therefore, defining the Debye relaxation time 



r D = (/2k B T 



(5.35) 



we have 



dg/dt = -g/r D 



g(t) = g(0) e-^ . 



(5.36) 



The initial condition g(0) = 1, comes from the matching of the distributions 
( I5.33P and (15.341) at t — 0. Since the normalisation constant follows from 
1 = Jq d$ sin $Po = 2N, we finally have 



(5.37) 



It is easy to compute now the average dipole moment along the field direction 
{p cos^) = d$ sm{}P({}, t)p cos by the change of variables z = cosfi 




(p cos 



V 



dz 



-i 



P AE /td2 

knT 




(5.38) 



This result goes from the Curie law for the linear response of a paramagnet 
p 2 AE/3k-BT in the initial equilibrium regime, to zero at long times, corre- 
sponding to the final equilibrium state in the absence of field. The solution 
of the Fokker-Planck equation (I5.37P provides also a complete description 
the intermediate non-equilibrium regime. 
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6 Linear response theory, dynamical suscep- 
tibilities, and relaxation times (Kramers' 
theory) 

In this section we shall consider some general results that can be obtained 
for the response of a very general class of dynamical systems (and in par- 
ticular those having a "Fokker-Planck dynamics"), results which hold when 
the external perturbation is weak enough. In the context of this linear re- 
sponse theory, due to Kubo, the definition of dynamical response functions 
appears naturally (in the time and frequency domain), and sometimes, asso- 
ciated with them, quantities characterising the lifetime of certain states — the 
relaxation times. 

6.1 Linear response theory 

Let us consider a system governed by an evolution equation of the type 



dtP = CP , 



(6.1) 



where £ is a linear operator (linear in its action on P), and P characterises the 
system. The natural example in this context is the Fokker-Planck operator 

£ = ~d y (a« •) + \dl (a™ •) (6.2) 

Well, let us apply an external perturbation, and separate the part of the 
Fokker-Planck operator accounting for the coupling with the perturbation, 
and the unperturbed part, which we assume to have a stationary solution 
P : 

£ — £,0 + £ext{t) , £ flPo — o . 

These conditions are quite general. For instance C could be the Liouville 
operator of a mechanical system, P$ the equilibrium Boltzmann distribu- 
tion, and £ cxt the part corresponding to the external potential. In the 
Fokker-Planck case £ cx t could be V^ t d v in the Klein-Kramers equation or 
l~ l dxiyl^-) in the Smoluchowski equation. However, the external perturba- 
tion does not need to be a force or a field. For instance, external modulations 
of the system parameters, like the bath temperature AT(t) are also possible; 
then in the Smoluchowski equation we will have £ ext oc ATd^. 
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If the perturbation is weak enough, we can write the deviation from the 
stationary state as P = P +p, and the evolution equation would lead to first 
order to 

dtPo+d t p = [£, + C cxt {t)} (Pq + p) ~ C , P +£ , p + C cxt {t)P , 
(we have disregarded C ext p). The resulting equation can be solved formall}@ 



d t p = C,oP + C cxt (t)P 



This equation gives the formal solution to the problem of time evolution in 
the presence of a weak time- dependent perturbation. ToDo, warn on the 
order of the operators 




(6.3) 



6.2 Response functions 
6.2.1 Time domain 

Let us consider any function c(y) of the variables of the system y, and com- 
pute the variation of its average with respect to the unperturbed state 



AC(t) = (c) (t) - (c) . 



(6.4) 



To this end we extract the time dependent part of C e xt(y, t) = C sx t(y)F(t) 
(factorisation is the common case) and use the solution (I6.3P 



AC(t) 



dyc{y)P{y,t)- ldyc(y)P (y) 
dyc(y)p(y,t) 



ds 



dyce^ C ' C ext P 



F(s) . 
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Indeed, taking the t derivative of the presumed solution, we have 



d tP = £cxtOO£b + C ,°J d S e {t - s)c - ( 'C cx t(s)Po=C cxt {t)P +C fi P Q.E.D. 



integrand at s = t 



£ ,op 
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Then, introducing the response function for the quantity c 



Rc(t) 



Jdyc(y)e tc - C cxt P (y) t>0 



(6.5) 



t < 



we can write the response AC(t) simply as 




(6.6) 



The following linear response functions are used [Q(t) is the step function 
e(t < 0) = and Q(t > 0) = 1]: 



S(t) pulse response function AC p (t) 
F(t) = ^ O(t) excitation function AC e (t) 
Q(—t) relaxation function AC r (t) 



(6.7) 



Let us consider the last one, also called after-effect function, which corre- 
sponds to switch a constant excitation off at t = 



/oo poo 
dsR c (t- s)e(-s) = / dsR c (t-s)= / ds' 
oo J —oo ^ s/ ^ Jt 



Rr(s') . 



Thus, we see that R c (t) can be obtained as the derivative of AC T (t) 



AC T (t) 



ds R c {s) 



R(t) = -—AC T . 
w dt 



(6.8) 



6.2.2 Frequency domain 

Introducing now the Fourier transforms of AC(t), F(t), and R c (t) 

AC(u) = fdte-' luJt AC(t), F(cu) = [ dt e~ iuJt F (t) , X c{u) = / dt e~ iwt R c {t) 



the convolution in Eq. f)6.6p relating those quantities, reduces to a simple 
product 



(6.9) 
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The complex quantity Xc{w) is known as the susceptibility. It can be seen 
that corresponds to the usual definition: if we excite with a perturbation 
F(t) = e luJt , the corresponding response function AC(t) oscillates in the 
stationary state with e luJt with proportionality coefficient x c (o;) 



AC(t) 



dsR c (t- s)e il 



ds R c {t — s)e 



-io>(t— s) 



Suit 



To conclude, we shall demonstrate a famous linear response relation be- 
tween the dynamic susceptibility and the relaxation (after-effect) function. 
From the definition of Xc{w), using that R c (t < 0) = and its relation with 
the relaxation function, we have 



Xc{UJ) 



def. 



dte~ 10Jt R c (t) 
-AC T (t)e~ iu}t 



dte 



—iuit 



-d(AC r )/dt 

RJf) 





oo 



+ / dtACJt)— e~' lwt 
o Jo dt 



AC r (0)-iw / dtAC v (t)e 
Jo 



-iuit 



Then, extracting the static (thermal-equilibrium) susceptibility x c c q — Xc(0) 
AC r (0), we have 



Xc(w) = X 



1U 



dt A^) 
AC r (0) 



(6.10) 



which gives the dynamical susceptibility in terms of the normalised relaxation 
function AC r (t)/AC r (0). 



Example: Debye relaxation. We can immediately apply these linear 
response results to the example of the dynamics of the electrical dipole [Eq. 
f l4.30p ]. There, we calculated the time evolution of the average of the field 
projection of the dipole, Eq. (15.381) . Thus, in this case 



c = p cos "d 
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Debye Relaxation 




0.001 0.01 0.1 1 10 100 1000 

wt 



Figure 6: Real part (full line) and imaginary part (dashed line) of the Debye 
formula (16. lip for the dynamic susceptibility (normalised). 



and hence Xp q = (p 2 /3k B T) (Curie law). Therefore, we have AP r (i)/AP r (0) 

e -t/TD go that 



1/(>+1/td) 



P 



3k B T 



l-ko dte~ (iw+1/TD) * 




which is the famous formula displaying "Debye" relaxation (see Fig. [H]). 



6.3 Relaxation times 

When we considered various examples of Fokker-Planck equations, we ob- 
tained the solution for the Smoluchowski equation of an harmonic oscillator 



P{x,t) 



27ck B T{l - e- 2t /\ 



exp 



x e 



-t/T\2 



2k B T{l - e- 2t / T ) 



T 



7/w, 



In this equation we see that the time scale for the relaxation to the equilib- 
rium state Po oc exp(— ^mu^x 2 / k B T) , is given by r = ■j/uJq. This quantity is 
the relaxation time. In this problem it depends on the system parameters 7 
and Uq, but it is independent of the temperature. 

Now, in the example of the dielectric dipole, a natural relaxation time 
has also appeared Tp = Q/2k B T [Eq. (15.351) ]. In this problem the relaxation 
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Figure 7: Sketch of a metastable potential for the Kramers calculation of the 
relaxation time. The point a is at the potential minimum, b at the maximum, 
and c is a point at the right-side of the barrier. 

time depends on T, which is the common case, however, the dependence is 
not very strong. 

It is very common to find expressions for different relaxation times that 
depend exponentially on T (Arrhenius law), which is the generic behaviour 
when to establish the equilibrium potential barriers need to be overcome. 
This was not the case of the previous examples, and such dependence was 
absent. We shall solve now a simple problem with potential barriers to see 
how the exponential dependence arises (the theoretical study of this problems 
was initiated by Kramers in 1940 to study the relaxation rate of chemical 
reactions). ToDo, warn on low temperature assumption 

To simplify the calculation let us consider an overdamped particle, de- 
scribed by Smoluchowski equation 



where the last equality defines the current of probability J, and the metastable 
potential is depicted in Fig. [7J At very low temperatures the probability of 
escape from the metastable minimum is very low (zero at T = 0; determinis- 
tic system). Therefore the flux of particles over the barrier is very slow, and 
we can solve the problem as if it were stationary. Then the expression for J, 




(6.12) 



- 7 J = V'P + 



k B T8P 

m dx ' 



(6.13) 
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is assumed to be independent of x and the differential equation for P can be 
integrated (U = mV, D = k B T '/Wy) 



-U/k B T 



c J 



dx ' e U^')/k B T 



(6.14) 



were c is an arbitrary point. The integration constant is C% = P(c)e u ^ c ^ kBT . 
If we choose c well outside the barrier region c — > oo, we have P(c) ~ 0, and 
we find for the current J 

j=- D _Q*E . ( 6 .i5) 

dx' e u ( x ')/ k v T 



Since dJ/dx ~ we can choose x at will; we set x = a (the metastable 
minimum), so the integral in the denominator covers the entire maximum. 
The main contribution to that integral comes from a small region about the 
maximum x = b, so we expand there U(x) = Ub — \mu 2 {x — b) 2 , being 
muJl = U"{b). The integration limits can be shifted to ±oo, so the resulting 
Gaussian integral leads 



J = D 



P( a ) e Ua/k B T 



e u b/kBT^ 2 7Tk B T/ 



DP(a) 



mUJ b ( Uh - Ua )/k B T 



mut 



27rk B T 



(6.16) 



To compute P(a) we use the following argument. The fraction of particles 
close to the potential minimum iV a can be obtained integrating P(x) in an 
interval around a, with the distribution approximated as P = C2e~ u ^ x ^ kBT 
with U(x) = U a + \mu\ix — a) 2 , where mu 2 = U"(> 



Then for the particles 
2 so that P(a)/N a = 



in the well we have iV a = C 2 e- Ua/k * T ^2nk B T/mu i 
\fmu^j2/Kk^T . 

The relaxation rate is defined as J/N a (so the number of particles in the 
well N a , times the escape rate, gives the flux J leaving the well). Then, 
introducing the above expression for P(a)/N a into Eq. (16.161) divided by N a , 
and using D / k B T = 1 / 7717 we finally have 



1 

T 



WgUb c -(U h -U a )/k B T 



27T7 



(6.17) 



This formula has the typical exponential dependence on the barrier height 
over the temperature. Although the calculation in other cases (intermediate 
to weak damping) is much more elaborated, this exponential dependence 
always appears. 
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7 Methods for solving Langevin and Fokker— 
Planck equations (mostly numerical) 

In general the Fokker-Planck or Langevin equations cannot be solved analyt- 
ically. In some cases one can use approximate methods, in others numerical 
methods are preferable. Here we shall discuss some of these methods. 



7.1 Solving Langevin equations by numerical integra- 
tion 

We shall start with methods to integrate the Langevin equations numeri- 
cally. These are the counterpart of the known methods for the deterministic 
differential equations. 



7.1.1 The Euler scheme 

In order to integrate the system of Langevin equations 

^ = ^(y,t)+]TB ifc (y,t)£ fc 0O , (£ fc (t)> = , = 2DS M 5(t-s) 

k 

starting at t = to with the values yo, to the time to + T, one first divides the 
time interval [t , to + T] into N s time steps of length At, i.e., t n = t + nAt. 
The stochastic variables at a later time y(t n+ i), are calculated in terms of 
y(t n ) according to 



J/i(*n+i) = Viitn) + af ] [y(t n ),t n }At + ^B ik [y(t n ) 1 t n }AW kn 



(7.1) 



where = A4 + D^ Jlt Bj k djB ik is the first jump moment, AW kn , k = 
1, . . . , A^l (the number of Langevin sources), n — 1, . . . , iV s , are independent 
Gaussian numbers with zero mean and variance 2D At, i.e., 

(AW kn ) = , {AW kn AW k , n> ) = (2DAt)5 w 5 nn , . (7.2) 

The recursive algorithm (17. ip is called the Euler scheme, in analogy with the 
Euler method to integrate deterministic differential equations. By construc- 
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tion, for At — > 0, the above recursive scheme, leads to the correct Kramers- 
Moyal coefficients!^! 



7.1.2 Stochastic Heun scheme 

This is a higher-order scheme for the numerical integration of the Langevin 
equations given by (a sort of Runge-Kutta scheme) 



Vi (t + At) = Vl (t) + §{A[y,* + At]+A[y(t),t]}A* 

+ | Yl i B *\y> * + M + B ik [y(t),t}} AW k , 

k 



with Euler-type supporting values, 

yt = Vi (t) +A i \y(t),t]At + ^2B ik \y(t),t]AW k . (7.5) 

k 

Note that if one uses this support value as the numerical integration algo- 
rithm [by identifying y^t + At) = y*], the result does not agree with the 
ordinary Euler scheme if dB ik /dyj ^ (or equivalently if op ^ Aj). 

The Euler scheme only requires the evaluation of Ai and Bi k at one point 
per time step, while the Heun scheme requires two, increasing the compu- 
tational effort. Nevertheless, the Heun scheme substitutes the derivatives of 
B ik by the evaluation of B ik at different points. Besides, it treats the deter- 
ministic part of the differential equations with a second-order accuracy in At, 

25 Let us prove this in the simple one-variable case. Then 

y(t + At) = y(t) + a« (y, t)At + B(y, t)AW . (7.3) 

To obtain the Kramers-Moyal coefficients, we average this equation for fixed initial values 
y(t) (conditional average). To do so, one can use (AW) = and (Alf 2 ) = 2D At, to get 

(BAW)=0, (a^AtBAW) = , (BAWBAW) = 2DB 2 At . 

Therefore, one obtains 

(y(t + At) - y(t)) = a^At 
( [y(t + At) ~ y{t)f ) = (a (1) ) 2 Ai 2 + 2a^AtB (AW) + B 2 (AW 2 ) = 2DB 2 At + 0[{At) 2 } , 

which lead to the Kramers-Moyal coefficients (|5.11|) via Eq. (|4.17[) . Q.E.D. 



(7.4) 
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being numerically more stable. Thus, the computational advantage of the 
Euler scheme, may disappear if it needs to be implemented with a smaller 
integration step (At) to avoid numerical instabilities. 

7.1.3 Gaussian random numbers 

The Gaussian random numbers required to simulate the variables AW, can 
be constructed from uniformly distributed ones by means of the Box-Muller 
algorithm (see, e.g., Ref. [31 p. 280]). This method is based on the following 
property: if T\ and r 2 are random numbers uniformly distributed in the 
interval (0, 1) (as those pseudo-random ones provided by a computer), the 
transformation 

W\ = a/— 2 ln(ri) cos(27rr 2 ) , w 2 = a/— 2 ln(ri) sin(27rr 2 ) , (7.6) 

outputs w\ and W2, which are Gaussian-distributed independent random 
numbers of zero mean and variance unity. Then, if one needs Gaussian 
numbers with variance a 2 , these are immediately obtained by multiplying 
the above Wi by cr (e.g., a = ylDAt in the Langevin equations). 

7.1.4 Example I: Brownian particle 

The Langevin equations for a particle subjected to fluctuations and dissipa- 
tion evolving in a potential V(x) are 

{dx I dt = v 
dv/dt = -V'--yv + Z(t), {Z(t)£(s)) = 2D6(t-s) (7 ' 7) 

with D = j(k-B,T/m). For the potential we consider that of a constant force 
field F plus a periodic substrate potential of the form 

V(x) = -d [sin x + (p/2) sin 2x] - Fx . (7.8) 

For the familiar case p = 0, we have the cosine potential and the Langevin 
equations describe a variety of systems: 
(i) Non-linear pendulum: 

+ 70 +(#/•<) sin0= N + £(t) . (7.9) 

In this case we have x = <fi (the angle of the pendulum with respect to the 
vertical direction), al = gji (gravity over length of the pendulum), F = N 
(external torque), and D = ^{k^T /m€). 
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(ii) Josephson junction (RCSJ model): 

CV + |0 + f /c sin = | [/ + £(*)] . (7.10) 

Here x = (p (the phase across the junction), 7 = 1/RC, d = 2eI c /hC (essen- 
tially the critical current), F = 2eI/hC (external current), and D = k^T / R. 

(iii) Others: superionic conductors, phase-locked loops, etc. 

When p 7^ in Eq. (17. 8p . V is called a ratchet potential, where it is more 
difficult to surmount the barrier to the left than to the right (like a saw tooth; 
see Fig. [8]). Ratchet potentials have been used to model directional motion 
in diverse systems, one of them the molecular motors in the cell. 

If Fig. Owe show the average velocity vs. force for a system of independent 
particles in a ratchet potential, obtained by numerical integration of the 
Langevin equation (17.71) with a fourth order algorithm. It is seen that the 
depinning (transition to a state with non-zero velocity), occurs at lower forces 
to the right that to the left, as can be expected from the form of the potential. 
It is also seen that for lower damping, the transition to the running state is 
quite sharp, and the curve quickly goes to the limit velocity curve 7 (v) = F. 
The reason is that for high damping, if the particle has crossed the barrier, 
it will not necessarily pass to a running state, but can be trapped in the next 
well, while the weakly damped particle has more chances to travel, at least 
a number of wells. 



V(x)=-d[sin(x)+0.22sin(2x]]-Fx, F-0,0.1 ,-0.1 

2 | 1 1 




-6 1 ' 1 ' 1 

-10 -5 5 10 

position 

Figure 8: Ratchet potential (I7.8p . with asymmetry parameter p = 0.44, at 
zero force and F/d = ±0.15 (displaced for clarity). 



68 




Figure 9: 7 (v) vs. F for a particle in a ratchet potential at T/d = 0.1 with 
7 = 0.5 and 5. The dotted line shows the limit curve 7 (v) = F. Results 
obtained by numerical integration of the Langevin equation with a Runge- 
Kutta-like 4th order algorithm for a system of 1000 independent particles. 



v 



5000 10000 15000 20000 

TIME 

Figure 10: Time evolution of the position of two independent Brownian 
particles in a ratchet potential corresponding to Fig. The particles are in 
a weak force F/d ~ 0.03 so their random- walk is biased in the force direction. 
The other parameters are T/d ~ 0.1, 7 = 0.5. 

The smooth graphs in Fig. [9] are obtained averaging the results for 1000 
particles. The individual trajectories of the particles, however, are quite 
irregular. In Fig. [TQ], the trajectories of two of them are shown. It is seen 
that to the overall trend of advancing in the direction of the force, there are 
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superimposed Brownian fluctuations (biased random walk), and indeed we 
see that the particles can be trapped in the wells for some time and even to 
return to the previous well. 

7.1.5 Example II: Brownian spins and dipoles. 

The Langevin equation for a spin subjected to fluctuations and dissipation 
is the Landau-Lifshitz equation 

-A?a(sAb) , (&(*)&(«)) = 2D6 kt 6{t-s) , (7.11) 



as 

— = s A 
dt 




Figure 11: Symbols: x(a;,T) vs. T obtained by numerical integration of the 
stochastic Landau-Lifshitz equation with the Heun scheme (17.41) . Thin lines: 
simple Debye approximation. Thick line: thermal-equilibrium susceptibility. 
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Figure 12: Projection of s(t) onto the anisotropy axis for 7i = —AUs 2 z in 
zero field, for various temperatures. 



with D = Xk B T/s???l [cf. D = j(k B T/m) for the particle] and B = -&H/ds 
(cf. —dV/dx). For the Hamiltonian we choose 

H(s) = -AUs 2 z - s ■ B , B = 2AUs z z + B , (7.12) 

with the term of coupling with the external field Bo and a magnetic anisotropy 
term (magnetocrystalline, magnetostatic, etc.), which favours orientation of 
the spin along ±z. The anisotropy and the field play the role of the substrate 
potential and the force for the particle problem. 

The stochastic Landau-Lifshitz equation describes a variety of systems 

(i) Magnetic systems: magnetic nanoparticles and approximately, mag- 
netic molecular clusters. 

(ii) Electric systems (taking formally the limit A ^> 1): permanent dipole 
molecules, nematic liquid crystals, and relaxor ferroelectrics. 

Figure [TT] displays the results for the dynamical susceptibility vs. the 
temperature for an ensemble of 1000 spins with parallel anisotropy axes for 
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A = 0.1. At low temperatures, the relaxation time obeys r ^> 2tt/uj, so that 
the probability of over-barrier rotations is almost zero. The spins rotate close 
to the bottom of the potential wells (see the panel ksT/AU = 0.12 of Fig. 
[T2"]) and the response is small. As T is increased the spins can depart from 
the minima (panel ksT/AU = 0.18) and rotate over the barrier. However, 
since the jump mechanism is not efficient enough, the response lags behind 
the probing field, and an appreciable imaginary component of the response 
arises. If T is further increased, the very thermal agitation, which up to this 
point was responsible for the growth of the response, reaches a level that 
simultaneously disturbs the alignment of the magnetic moments in the field 
direction. The response has then a maximum and starts to decrease. Finally, 
at still higher temperatures the jumps are so frequent (panel k B T/AU = 0.4) 
that the spins quickly redistribute according to the conditions set by the 
instantaneous field and the response tends to the thermal equilibrium one. 

7.2 Methods for the Fokker— Planck equation 

There are several methods to solve Fokker-Planck equations not amenable 
for an analytical treatment. Among them we can mention the method of the 
eigenvalues, which consist in finding the eigenvalues of the Fokker-Planck 
operator occurring in the dynamical equation d t P = C FP P However, the 
operator C FP is not in general Hermitian. If the stationary distribution is 
known P , and L FP fulfills the detailed balance condition, it can be seen that 
the transformed operator C FP = P^ 1 ^ 2 C FP P^ 2 is Hermitian, so the problem 
can be reduced to an ordinary eigenvalue problem. 

£fpPu = ~KPn , (7.13) 

from which we will have P(y, t) = a n p n (y)e~ Xnt . However, along with re- 
lying on the knowledge of Pq and the detailed balance condition, the method 
has the drawback that the eigenvalue problem may be difficult to solve (as 
happens for the Schrodinger equation). Concerning approximate methods, 
we can mention the small noise expansion (also for Langevin equations) and 
the method of adiabatic elimination of fast variables, which reduces the di- 
mensionality of the problem. 

In what follows we shall discuss an alternative method, the continued 
fraction method, which is a special case of the expansion into complete sets 
approach. The method, although limited in principle to problems with a few 
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variables, it is exact and in general very efficient, and, in addition, illustrates 
more general techniques used in non-equilibrium statistical mechanics. We 
shall discuss the method with two important examples: the Klein-Kramers 
equation and the Fokker-Planck equation for the Debye dipole. 

7.2.1 Solution of the Klein— Kramers equation by continued— fractions 

The Klein-Kramers equation (15.201) can be written in the compact form as 

/ \ I A-ev = — v d x + V'd v = {H , •)__ 

d t P= ( C ICV + C irr )P , I 1 /PB (7.14) 

Scaled quantities. To simplify the notation, we introduce a thermal rescal- 
ing of the velocity, time, damping, and potential 

v = v/ ^/k B T/m , i = t x a/ k B T/m , j — j/y/ k B T/m , U = mV/k B T. 

(7-15) 

Well, for the v, t, and 7 so defined we shall not keep the bars and we shall 
simply write the Klein-Kramers equation as 



d t P — (C T ev + Arr) P 



£ rev = -v d x + U'd v 

Arr = jd v (v + 8 v ) 



(7.16) 



Expansion in an orthonormal basis of v. We can expand the proba- 
bility distribution P(x,v,t) in an orthonormal basis ip n {v) (to be specified 
below) as follows 

P = Wj2cn(x,t)Mv) , (7.17) 

n 

where W = W(x,v) is some function we extract for later convenience. Due 
to the orthonormality of the if) n {v), we can write c n = Jdvip n P/W, and 
hence 

OBS: partial derivatives^ = J ' Av^W' 1 = J dv ip n {W~ l C FP W) PJW^ 
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Then, we can compactly write 



Cn ^ Qn,m,C"iri j 


Qn,m = j dv 1p n C FP ll) m C FP = W 1 C FP W . 


m 





(7.18) 

Then, the equations for the c n are like a system of linear equations with the 
matrix elements of £ FP . However, the Q n , m still contain operators over x. 
Besides, it seems that all the c m , contribute to the equation for c n . 

As for the function W, by analogy with the transformation P C FP P 
to simplify the problem, we try it at least for the v dependent part. Then, 
since Pq oc exp(— v 2 /2) in our unit@ we put W cx exp(— 1> 2 /4). For the x 
dependent part, one sets W oc exp(— ell), which for e = 1/2 corresponds to 

1 /2 

the x-dependent part of P (if given by the Boltzmann distribution). Thus, 
from these considerations we set W oc exp[— (v 2 /4 + ell)}. 

Calculation of C FP = W~ x C FP W . Since £ FP = £ rcv + C IVV) we have to 
calculate £ rev = W' x L reY W and £ irr = W~ x Ci„W. 
For £ irr we have 

7 -1 Arr/ = e*'% (v + d v ) e"^ 4 / = - {-dl + \v 2 - |) / , 

which has the structure of (minus) the Hamilton operator of the harmonic 
oscillator in quantum mechanics. Hence, we introduce creation and annihi- 
lation operators b + and b, defined as 

{I'-itt [M+] = 1 ' (7 ' 19) 

Since the 1/2 in £ irr cancels the "zero-point fluctuations", we finally have 

£ irr = _ 7 6+fe . (7.20) 



Because of this, it seems natural to choose the if) n (v) as the Hermite functions 

Mv) = r^—-r 2 e- v2/i H n (v) . (7.21) 



Undoing the thermal rescaling (|7.15|) . we have exp(— v 2 /2) — > exp(— mv 2 /2kv,T). 
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For £ rev , since £ rcv = — v d x + U'd v , we need d v and d x 



d x = d x - £[/' . 



= e £C/ ^ (e- £l/ /) , 
Then, for £ rev we have 

Cm = - {b + b + ) (d x - eU'l -b+U' = -b (d x - eU') - b + [d x + (ls)U'] 

v d v 

Calling D + to the "coefficient" of b and LL to that of b + , we finally have 



D+ = d x -eU' 
D_ = d x + (l-e)U' 



C rcv = -bD + - b + D_ . 



(7.22) 



Calculation of Q n ,m- the Brinkman hierarchy. Recall that we need 
£ FP = £ rev + £i rr , to obtain its matrix elements <3n,m = J dvip n jC FP ip m , which 
enter in the equation of motion for the coefficients c n = J2 m Qn,mC m - This is 
an easy task now, since we have written £ irr and £ rev in terms of b and b + . 
Then, using the "ladder" properties of b and b + , that is, b + ip n = \/n + lip n +i 
and hpn = \fnip n -i, plus the orthogonality of the ip n , we obtain 

Qn, m = Jdvtpn (£ rcv + £ irr ) = -VnD_S n -i !Tn -njS n>m -y/n + lD + 5 n+ltTn 

Therefore, the sum in c n = J2 m Qn,mC m , is trivially done, and one finally gets 
the so-called Brinkman hierarchy (1956) 



c n = - (y/nD- c n -i + 7nc n + Vn + 1D + c n+1 ^j . 



(7.23) 



We see that, due to (i) the choice of the basis functions to expand P and 
(ii) the extraction of the factor exp(— v 2 /4), only the nearest neighbours of 
c n contribute in c n = J2mQn,mC m - Writing explicitly the equations, we see 
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that this fact results in a tri-diagonal structure of the system 



-co = 


07C0 


+ 




+ 


+ 








-Cl = 




+ 


17C1 


+ V2D + c 2 


+ 





+ 




-c 2 = 





+ 


V2LLci 


+ 2 7 c 2 


+ 


V3D + c 3 


+ 





-C3 = 




+ 







+ 


37C3 


+ 





This equation is completely equivalent to the Klein-Kramers equation ( 17.161) 
and valid for any potential. 

Continued fractions. Why are we so happy for having transformed the 
Klein-Kramers equation into an infinity hierarchy for the c n ? First, by taking 
the Laplace transform /(s) = J °° dt e~ st f(t) a differential- recurrence relation 
of the general form 

dc- 

+ Q-a-x + Q iCl + Qfc i+1 = fi , (7.24) 
can be reduced to (we omit the tildes on the Laplace transforms) 

QiCi-t + QiCi + Qfc i+1 = fi , (7.25) 



where Qi = Qi + s and = + Cj(0). Then, this relation can be solved by 
introducing the ansatz c, = SiQ-i + a i: obtaining 



C -i-U c Qi Qt a i+i ~ fi 

Ci = + <2j with bi - 



Qi + Qf^i+l Qi + Qf Si+\ 



(7.26) 

It is to be remarked that the quantities involved in the relation (I7.25P do not 
need to be scalars, but they can be vectors (q and fi) and the coefficients 
Qi matrices. The only change in the solution (I7.26P is that the fraction bar 
then means matrix inversion. 
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The reason for the name "continued fraction", is that, if we consider for 
instance Si, it is given in terms of S^ + i in the denominator. But this can be 
expressed in terms of Si+2, and so on, leading to an expression of the form 



K 



Pi 



(7.27) 



9i + 



92 + 



93 H 



which is called a continued fraction. Well, if the q are judiciously chosen to 
decrease with increasing i, we can truncate at some large N , setting = 0. 
This leads all the quantities to vanish at i = N, and we can iterate downward 
the continued fractions in Eq. (I7.26P down to % = 0, storing the successive Si 
and aj. Then, starting from cq (usually known by some means, e.g., normali- 
sation of the distribution), we iterate upwards with Cj = SiCi-x+cii, obtaining 
the solution to the recurrence-relation (17.251) . To ensure convergence, one can 
repeat the calculation with a truncation index 2N, AN, etc. and check that 
the results do not change. Note that the steps described above are very easy 
to code in a computer program (even in a programmable pocket calculator). 
This is the answer to the question of why it was important to transform the 
Klein-Kramers equation into the Brinkman hierarchy. 

Solving the Brinkman hierarchy. As written, the hierarchy (17.231) . in- 
volves coefficients that still contain the operators D± over the x variables, 
instead of true coefficients. But we only need to find a representation of those 
operators in a basis of functions of ). Then we just calculate the 

matrix elements of any operator A in that basis A pq = J dxu*Au q , while the 
expansion coefficients of c n (x, t) are expressed as a column vector 



Then, the Brinkman hierarchy f 1 7 . 2 3 1) is reduced to the following differential 
recurrence relation 





C 



n 



P 



\ < J 



c 



n 



Q n C n _i + Q rt C n + Q„ C n +i 



(7.28) 
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whose coefficients are not any more operators but ordinary matrices, with 
elements 



(Qn)pg 
(On), 

(Q^ 



q 
pq 



'n 

-njSpg 



p'i 



(7.29) 



-v^Tl {d x ) pq -eU' pi 



The choice of the basis functions u p (x) is naturally dictated by the symmetry 
of the potential. 

In this form, we can directly use the method of continued fractions to solve 
the Brinkman hierarchy. In principle, with the c v n obtained we can construct 
the probability distribution and compute any observable. Nevertheless, this 
is not even required, since common observables are directly obtained from 
the expansion coefficients. For instance, if e = in W oc exp(— v 2 /A), this 
quantity is after normalisation W — ipo (see the definition of the Hermite 
functions (17.211) ] so that for the averaged velocity (v) = Jdx J dvvP(x,v,t) 
we get 



( v ) = J2 c ™ dxu p( x ) / 

np 



ipl 

dv ^_v^ipo{v) ip n (v) 

b+b+ 



J dx u p (x) 



where he have just assumed that the Oth element of the basis u p (x) is constant 
(say Uq = 1), so that the result follows from orthogonality (w -L u n ). 



Example: particles in a periodic potential. Let us consider the case 
of a general periodic potential U' = U' q e iqx , where the U' q are the Fourier 
components of the potential derivative. In this case the natural choice of 
the basis functions u p (x) is that of plane waves u p (x) = e igx /V2ir. Then 
the matrix elements needed to construct the matrices Q n in Eq. (I7.29P are 
simply 

@.) P9 = &pi, (U') pq = U p _ q . (7.30) 
Then, the number of diagonals occupied in the matrices Q n below and above 
the main diagonal is equal to the number of harmonics in the potential. In 
the example of the cosine potential, this number is only one, while in the 
ratchet potential (17.81) . this number is two. 

Well, in the example of the ratchet potential we computed the average 
velocity vs. force by Langevin simulation (Fig. [H]). If we compute now the 
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Figure 13: Depinning curves in a ratchet potential of Fig. [9j obtained with 
Langevin simulation (symbols), together with the results obtained with the 
continued-fraction method (thick lines). 

same curves with the continued-fraction method, we get a complete agree- 
ment shown in Fig. (T3] (without any free parameter). 

As an additional example, in Fig.[14l we show the real part of the dynam- 
ical susceptibility as a function of the frequency in an ordinary (non-ratchet) 
periodic potential (thus, the curves could correspond to a Josephson junc- 
tion, or to the ac conductivity in a transport problem). We see a peak about 
the frequency of oscillations in the bottom of the potential wells (which de- 
fines the frequency scale in the graph). The peak becomes more and more 
sharp and high the lower the damping is, showing the typical features of a 
resonant system. 

When the continued-fraction method can be used its advantages are: (i) 
it has no statistical error bars, (ii) it is in general extremely efficient, (iii) it 
outputs the complete distribution if required, and (iv) it may be extended 
to problems of quantum Brownian motion. The drawbacks are that it is 
quite specific of the problem to be solved (we need to compute the matrix 
elements for each potential, and the relation between the coefficients and 
the observables) , and in this respect the Langevin simulations do not have 
this problem. Besides, the simulations output trajectories, which are lost 
in the continued fraction approach. Finally, it can only be implemented for 
systems with a few variables (e.g., independent particles), while the Langevin 
simulations do not suffer this limitation. 
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resonance in ac field 
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Figure 14: Real part of the dynamical susceptibility (ac conductivity) vs. fre- 
quency. The curves have been obtained with the continued-fraction method 
for T/d = 0.1 and various values of the damping 7 = 1, 0.3, 0.1, 0.03. The 
resonance at u ~ 1 is clearly seen. 



7.2.2 Solution of the Debye equation by continued— fractions 

This continued-fraction method can also be applied to rotational problems. 
We shall illustrate this with the Fokker-Planck equation for a dipole f)4.30p . 
First, using the definition of the Debye relaxation time = (/2k-gT [Eq. 
f !5.35p ]. the abbreviation a = pE/k^T and changing to Cartesian coordinates 
z = cos$, the Debye equation can be written as 



2r dP _ d 

dt dz 




(7.31) 



In this problem the natural choice for basis function is the Legendre 
polynomials p n ( cos $) [cf. Eq. (17.171) 



P = ^2c n (t)p n (z) . 



(7.32) 



Due to the orthogonality of the p n (z), f_ x dz p n (z)p m (z) = 2<5 n m /(2n + 1), 
we have c n = \{2n + l)/2] J dzp n (z) P(z). Then the equation of motion for 
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□ebye relaxation in a dc field 
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Figure 15: Dynamical susceptibility (polar isability) of a dipole in various 
static fields. The curves have been obtained with the continued-fraction 
method for T = 0.1 and the fields used are pE/k^T = 0, 1, 2 and 4. 



the c n reads [cf. Eq. fTT^gj) 

dtP 



4r D 
2n + 



^C-n — j dz p n -£-fp — ^ ^ Qn,m^rn i with Qn,m — J 



dz p n £, FP p ri 



To calculate Q njm (which is not an operator here), we use relations obeyed 
by the Legendre polynomials^! and an integration by parts: 



Q, 



dzp r , 



dz 



(1 - z 2 ) 



dp m 
dz 



~ ap n 



-m(m + 1) J dzp n p m + a Jdz - z 2 )p, 

n(n + 1) / 25 n _i >m 



-n(n + 1 



25 nm 

2n + l 



+ a 



2n + l \2n-l 



27 



Specifically, the Legendre equation, and a sort of first integral: 



_d 

dz 



d^ 

dz 



on d»„ n(n + 1) 
+ n(n+l)p n = 0, (l-^)^L = ^_^_i 



2n + 3 



(Pn-l -Pn+l) 
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Inserting this result into the equation for c n , we finally find 



2r D 




pE ( c n _i 



c n+l 



) 



(7.33) 



n(n + 1) 



k B T \2n - 1 



2n + 3 



which is a tridiagonal differential-recurrence relation for the dipole analogous 
to the Brinkman hierarchy (17.231) . 

Solving Eq. (17.331) by continued fractions, we can obtain any observable. 
For instance, we can compute the linear susceptibility in a static field (re- 
member that the solved this problem in the absence of bias field) , by setting 
E = E + AEcos(u)t). The results are shown in Fig. [151 f° r various Eq. 
It is seen that the equilibrium susceptibility (real part at uj — > 0) decreases 
with E , since the susceptibility measures the slope of the static polarization, 
which saturates at large fields. Besides the relaxation time (roughly, the in- 
verse of the location of the maximum of the imaginary part) decreases with 
Eq, so that the relaxation is faster in a static field. 
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8 Derivation of Langevin equations in 
the bath-of-oscillators formalism 



Now we shall discuss the foundations of the Langevin equations by showing 
how they can be "derived" from more or less microscopic description of a 
system coupled with its environment. 

8.1 Dynamical approaches to the Langevin equations 

The Langevin equations we have discussed for Brownian particles and spins, 
are phenomenological inasmuch as they constructed by augmenting known 
phenomenological equations by fluctuating terms. For subsequent reference, 
let us first rewrite these equations in the notation we shall employ later: 

• Brownian particle 

{dq p 
dt m 
it = -f - 1 Tt + ««> • <«*»« = WMW - • 

(8.1) 

Note that the fluctuating terms only enter in the equation for the momentum. 

• Brownian spins and dipoles 

- X s A (s A B) , (&(t)6(t')> = 2(Xk B T)5 M 5(t - t') . 

(8.2) 

This stochastic Landau-Lifshitz equation is equivalent to the Gilbert equation 
where the damping term is oc —As A (ds/dt). Recall that Eq. (18.21) contains 
as a special case the stochastic dynamical equation for the electrical dipole. 

Note that in both equations the fluctuating and dissipation terms are not 
independent: D = 7/c B T/m for the particle, which corresponds to D — Xk B T 
for the spin. Besides, the force F = —dV/dq, corresponds to B = —&H/ds, 
while the damping — j(dq/dt), corresponds to —As A ^s A Bj , which as men- 
tioned above is related with oc —As A (ds/dt). 

There have been several attempts to justify, starting from dynamical de- 
scriptions of a system coupled to its surroundings, these important Langevin 
equations. The effort was first directed to the Langevin equations for the 



as 



s A 



S3 



Brownian particle (translational problems) and then to rotational Brownian 
motion of spins and dipoles. 

In most of those studies, the environment is represented as a set of inde- 
pendent harmonic oscillators. The oscillators are somehow "projected out" 
and an equation for the system variables is derived. The final equation 
has the form of a generalized Langevin equation (i.e., containing "memory" 
terms), whose fluctuating and dissipative terms naturally obey fluctuation- 
dissipation relations. For instance, for the particle problem one gets 



dp m 

dt " ~ ~dq 



+ /(f) - J\t' /C(f - 0^(0 , (/(*)/(0> = **TK(t - t') . 

(8.3) 

Thus the relaxation (damping) term, involves a memory integral taken along 
the past history of the system. Besides, the memory kernel JC(t — t') is 
determined by the correlation properties of the fluctuating force /(f). Thus, 
if the autocorrelation of the fluctuating terms is very short (f(t)f(t')) oc 
5(t — t'), the damping term reduces to minus the velocity —(dq/dt) of the 
particle. Similar results are obtained for the spin and dipole problems. 

In what follows we shall discuss the bath of oscillators formalism. First, 
since we shall use a Hamiltonian formalism throughout, we shall start with 
a brief review of the main results from Hamiltonian mechanics that we shall 
need. Subsequently, we shall introduce the model for the system coupled 
with its environment, deduce the corresponding dynamical equations, and 
finally discuss some examples. 

8.2 Quick review of Hamiltonian dynamics 

The dynamical equations for a system with Hamiltonian T~C(p, q) are 

dg &H dp dH 



dt dp ' dt dq 



(8.4) 



Then, the time evolution of an arbitrary dynamical variable of the system 
A(p, q) (assumed not explicitly time dependent), is dA/dt = (dA/dq)(dq/dt) + 
(dA/dp)(dp/dt). Then, using for dq/dt and dp/dt the Hamilton equations 
and introducing the Poisson bracket of two arbitrary dynamical variables 

I A £l = — — -— — (8 5) 

^ ' -* dq dp dp dq 
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we have for dA/dt the basic Hamiltonian evolution equation 




(8.6) 



Finally, by using this equation for A = q, with dq/dq = 1, dq/dp = 0, 
and for A = p with dp/dq = 0, and dp/ dp = 1, we see that the Hamilton 
equations (18.41) are a particular case of Eq. (18. 6p . For a system with variables 
(.Pa, la) a = 1, • • • , N, the above results are the same, with the only change 
of introducing a sum over a in the definition of the Poisson bracket. 

Two more results we shall need are the product rule of the Poisson bracket 



{A,BC} = {A,B}C + B{A,C} , 



5.7) 



and the chain rule 



i . k 



Xi(p,q) 



which immediately follow from the ordinary differentiation rules. 



Spin dynamics case. The equations for an isolated classical spin (not 
subjected to fluctuations and dissipation) 



— = sAB 
dt 



B 



&H 

ds 



(8.9) 



can also be written in Hamiltonian form. To this end, let us write the formula 
for the gradient operator in spherical coordinates 

du ^du * 1 du A 1 du 
-t^ = s— + ¥- 

OS 



ds 



s sin dcp 



[8.10) 



where tp and $ are the azimuthal and polar angles of s. Since the length of 
s is constant, the set vectorial equations (18.91) . can be written as 



dip 
~dt 



1 dU 

s sin •& dd 



dd 



1 dU 



dt s sin $ dip 



(8.11) 



which correspond to the Hamilton equations (18 ,4p with the conjugate canon- 
ical variables 

q = ip , p = s z . I (8.12) 
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In terms of the variables (I8.12p the Cartesian components of the spin are 
given by 



s x = \A 2 — P 2 cos Q 5 s y = V s2 ~ P 2 sin g , s z = p . 

From these Si(p, q) and the definition of the Poisson bracket of two arbitrary 
dynamical variables [Eq. (18.51) ]. one can readily obtain the customary Poisson- 
bracket ( "commutation" ) relations among the spin variables 

{Si, Sj} =^ tijkSk , (8.13) 
k 

where e^fc is the Levi-Civita symbol@ In addition, on using the chain rule 
of the Poisson bracket [Eq. (18.81) ]. one gets the useful relation 

{ St ,W(s)} = -(sA^r) , (8.14) 



ds 

which is valid for any function of the spin variables W(s)^ 2 '' 



8.3 Dynamical equations in the bath-of-oscillators for- 
malism 

We shall now study a classical system surrounded by an environment that can 
be represented by a set of independent classical harmonic oscillators. In spite 
of its academic appearance, those oscillators could correspond to the normal 
modes of an electromagnetic field, the lattice vibrations of a crystal (in the 
harmonic approximation), or they can be an effective low-energy description 
of a more general surrounding medium (Caldeira and Leggett, [T]). 

28 To illustrate, from 

ds x /dq = — [s 2 — p 2 ] ^ 2 sing , ds x /dp = — p [s 2 — p 2 ] ^cosg, 
dsy/dq = [s 2 — p 2 ] ^ 2 cosq , dsy/dp = — p[s 2 — p 2 ] 1 ^ 2 sing, 

one gets {s x ,s y } = psin 2 q +pcos 2 q = s z . Q.E.D. 

29 Note that one can conversely postulate the relations {st,Sj} = ^2^^ijkSk, and then 
derive Eq. f|8. 9|) starting from the basic Hamiltonian evolution equation dsi/dt = {si,Ti\ 
and using Eq. (|8. 14|) . This can be considered as a justification of the presence of the 
expression B = — dH/ds in the dynamical equations for a classical spin. 
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8.3.1 The system-environment Hamiltonian 



The total system consisting of the "system of interest" plus the oscillators 
representing the environment forms a closed dynamical system that we shall 
describe by augmenting the isolated-system Hamiltonian as follows 

n T = n(p, g) + \{ p « + «£ [q° + ^«(p. <?)] '} ■ ( 8 - 15 ) 



a 



Here, a is an oscillator index [e.g., the pair (k, s) formed by the wave- vector 
and branch index of a normal mode of the environment], and the coupling 
terms F a (p,q) are arbitrary functions of the system variables. These terms 
may depend on the parameters of the oscillators ui a , but not on their dy- 
namical variables P a , Q a . On the other hand, we have introduced a system- 
environment coupling constant e for the sake of convenience in keeping track 
of the orders of the various contributions. 

The terms proportional to Fj, which emerge when squaring Q a +(s/uj^ t )F a , 
are "counter-terms" introduced to balance the coupling-induced renormaliza- 
tion of the Hamiltonian of the system. The formalism takes as previously 
considered whether such a renormalization actually occurs for a given inter- 
action, so that 7i would already include it (whenever exists). An advantage 
of this convention is that one deals with the experimentally accessible en- 
ergy of the system, instead of the "bare" one, which might be difficult to 
determine. 



8.3.2 Dynamical equations 

Let us first cast the Hamiltonian (18.151) into the form 

H T = wWfo q ) + J2 \ {Pi + u*Ql) + e Q«F a {v, q) , (8.16) 

a a 

where q and p are the canonical coordinate and conjugate momentum of a 
system with Hamiltonian T~C(p, q) and the "modified" system Hamiltonian 
augments Ji by the aforementioned counter-terms 

" a a 

Besides, in the above expression for 7ix the Hamiltonian of the oscillators is 
clearly recognised He = J2 a \ ( P a + w aQa)i an d the same for the coupling 
term H mt = e J2 a Q a F a (p, q)- 
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The equation of motion for any dynamical variable C without explicit 
dependence on the time, dC/dt = 0, is given by the basic Hamiltonian 
evolution equation (I8.6P with TC — ► TCt, with the whole Poisson bracket is 
given by 

, , _ dAdB dAdB ^ dA dB dA dB 
' J dq dp dp dq dQ a dP a dP a dQ a 

Therefore, the (coupled) equations of motion for any dynamical variable of 
the system A(p, q) and the environment variables read (C = A, P a , and Q a ) 

^ = {A,H^}+eJ2Q«{ A > F *} > ( 8 - 18 ) 

a 

^ = P a , ^ = -uZQa - eF a , (8.19) 

where we have used {Q a , Pa} = 1- The goal is to derive a dynamical equation 
for A(p, q) involving the system variables only {reduced dynamical equation). 

On considering that in Eqs. (18.191) the term — sF a (t) = —eF a [p(t),q(t)] 
plays the role of a time- dependent forcing on the oscillators, those equations 
can be explicitly integrated, yielding 

Qait) = Q h a (t) - — fdt' sm[u a (t - t')]F a (t') , (8.20) 

where 

Qa(t) = Qa(t ) COs[uJa(t - t )] + [P a (t )/uJ a } S in[uJ a (t - t )] , (8.21) 

are the solutions of the homogeneous system of equations for the oscillators 
in the absence of the system-environment interaction (proper modes of the 
environment). Then, on integrating by parts in Eq. (I8.20p one gets for the 
combination sQ a that appears in Eq. (18. 18j) 

(tf/C a (t- 0-^(0, (8-22) 



where 



f a (t) = eQl(t) , £ a (r) = 4 cos(u a r) . (8.23) 
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Next, in order to eliminate the environment variables from the equation 
for A(p, q), one substitutes Eq. (I8.22p back into Eq. (18.181) . getting 



dA 
~dt 



{A,H {m) } ^{A,F a }/C a (0)F a (t)+^{AF a }X: a (t-to)Fa(*o 



dt'K, a {t-t')^{t') 



The term ^ a [A, F a }K, a (t — t )F a (t ) depends on the initial state of the 
system (p(t ), q(t )) and produces a transient response that can be ignored 
in the long-time dynamics (we shall return to this question below). The 
parallel term — ^ Q {A, F Q } /C a (0) F a {t) is derivable from a Hamiltonian and 
balances exactly the term due to the counter-terms in {A, TC^'j. This can 
be shown by using 

- {A F a }JC a (0)F a = {A, -\ /C Q (0)F a 2 } , 



which follows from the product rule (18. 7p of the Poisson bracket and then 
using /C Q (0) = £ 2 / w q [ see Eq. (I8.23P ]. Therefore, one is left with the reduced 
dynamical equation 



^ = {An} + J2 (A F a } \f a (t) + fdt' K a {t - 1 
m n L Jto 



i\ dF a , , 



dt 



[8.24) 



where the first term yields the free (conservative) time evolution of the sys- 
tem, whereas the second term incorporates the effects of the interaction of 
the system with its environment. 

To conclude, let us decompose the coupling functions as 



[8.25) 



Here "a" stands for a general index depending on the type of interaction. 
The idea is to split the part of the coupling W a (p, q) which is common to all 
the modes, so that F a is obtained multiplying that part by certain mode- 
dependent system-environment coupling constants c^. For instance, if a is 
a mode of wave-vector k, and Ft = k ■ r , then a = i (the Cartesian index) 
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with c 1 -. = ki and Wi — rj. Introducing the above coupling function into Eq. 
( K2M . we have 



^ = {A,n} + J2 {A £ c* a W a } \f a (t) + fdt> K a {t - t>) <Z^W 



dt 



Therefore, we finally have 



^ = {a,h} + J2 (A w a } \f a (t) + fdt' V /c ab (t - 1 



dt 1 J 



where 




(8.26) 



(8.27) 



The terms f a (t) are customarily interpreted as fluctuating "forces" (or "fields") 
Indeed / a (£) is a sum of a large number of sinusoidal terms with different fre- 
quencies and phases; this can give to f & {t) the form of a highly irregular 
function of t that is expected for a fluctuating term (see below) |fj The inte- 
gral term keeps in general memory of the previous history of the system, and 
provides the relaxation due to the interaction with the surrounding medium. 

The origin of both types of terms can be traced back as follows. Recall 
that in Eq. (I8.20j) the time evolution of the oscillators has formally been 
written as if they were driven by (time-dependent) forces —eF a [p(t'),q(t')} 
depending on the state of the system. Therefore, Q a (t) consists of the sum 
of the proper (free) mode Q\{t) and the driven- type term, which naturally 
depends on the "forcing" (state of the system) at previous times. Then, the 
replacement of Q a in the equation for the system variables by the driven- 
oscillator solution incorporates: 

1. The time-dependent modulation due to the proper modes of the envi- 
ronment. 



30 Explicit expressions for the / a and the kernels in terms of the proper modes are 

Mt)=e^2<Qa(t) , U(r)^ 2 ^^cos( Wa r) (8.28) 



90 



2. The "back- reaction" on the system of its preceding action on the sur- 
rounding medium. 

Thus, the formalism leads to a description in terms of a reduced number of 
dynamical variables at the expense of both explicitly time-dependent (fluc- 
tuating) terms and history-dependent (relaxation) terms. 

8.3.3 Statistical properties of the fluctuating terms 

In order to determine the statistical properties of the fluctuating sources 
f a (t), one usually assumes that the environment was in thermodynamical 
equilibrium at the initial time (recall that no statistical assumption has been 
explicitly introduced until this point): 

Peq(P(*o), Q(*o)) oc exp [-/3W E (*o)] , n E (t ) = Yl \ [ p «(*°) 2 + "iQoihf] ■ 

a 

The initial distribution is therefore Gaussian and one has for the first two 
moments 

(Qa(t )) = , (Pa(t )) = , 

(Qa(t )Qp(to)) = 5 a p\- , (Q a (t )Pfs(t )) = , (P a (t )Pp(t )) = 5 a pk B T . 
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From these results one readily gets the averages of the proper modes over 
initial states of the environment (ensemble averages): 



(Qa(t)) = (Qa(to)) COsMt - t )] + CP«(*o)> ~ Sm[uJ a {t - t )] , 


(Ql(t)Q h a{t')) = (Q a {to)Qp{to)}cos[u a {t-t ))cos[u^t' -t )) 
v ^ , 

+ (Qa(to)P(3(to)) cos[u a (t - t )] sin[up(t' - t )] 
" v 'm b ujp 



+ {P a {t Q )Q p {t Q )) — sin[u; Q (t - t )] cos^t' - t )] 



+ (P a (to)Pp(to)) 



1 



■ u a m b U(3 



sm[u a (t - t )] sm[ujp(t' - t )] 



k B T ^{ cos [u a (t - t )] cos[^ Q (t' - t )] 



so that 



+ sin[a; a (t - t )] sm[u a (t' - t )]} 



(Q h a (t)) = , (Q h a (t)Q h ,(t')) = cos[u; Q (t - t')] . (8.29) 

Then, since Eq. (JH^HD says that / a (t) = ££ a c£Q*(t) and /C ab (r) = 
^Ea^^a/^a) 008 ^ 1 ")! the equations (18. 291) give for the averages of the 
fluctuating terms / a (t): 



(/.(*)> 

(/.(*)/b(0> 



o 



k B TJC ah (t-t') . 



(8.30) 



The second equation relates the statistical time correlation of the fluctuat- 
ing terms / a (t) with the relaxation memory kernels /C a b(r) occurring in the 
dynamical equations {fluctuation- dissipation relations). Short (long) corre- 
lation times of the fluctuating terms entail short-range (long-range) memory 
effects in the relaxation term, and vice versa. The emergence of this type of 
relations is not surprising in this context, since fluctuations and relaxation 
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arise as different manifestations of the same interaction of the system with 
the surrounding mediumjlj] To conclude, we show in Fig. [TH], the quantity 
/(O = £ ^2k c k [Qk(to) cos(u k t) + [P k (t )/uj k ] sin(u k t)), with c k oc k, u k = ck, 
and the (P(t ), Q(t )) drawn from a Gaussian distribution. The graph shows 
that a quantity obtained by adding many sinusoidal terms with different 
frequencies and phases can actually be a highly irregular function of t. 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

TIME 



Figure 16: The quantity f(t) obtained by summing over 1000 oscillators 
with initial conditions drawn from a Gaussian distribution. 

8.3.4 Markovian regime 

We shall now study the form that the dynamical equations derived exhibit 
in the absence of memory effects. This occurs when the memory kernels 

31 If one assumes that the environment is at t — to in thermal equilibrium in the presence 
of the system, which is however taken as fastened in its initial state, the corresponding 
initial distribution of the environment variables is P eq oc exp[— Hse^oV^bT 1 ], with 

H SE (t Q ) = J2 a H p M 2 +u; 2 a [Q a (to) + (e/ujDFM] 2 } . 

In this case, the dropped terms JC a (t — t )F a (t ), which for F a = ^2 a c%W a lead to 
J2b fc&b(t— to)Wh(to), are included into an alternative definition of the fluctuating sources, 
namely f a (t) = f a (t) + J2h ^ab(* — to)Wb(to). The statistical properties of these terms, as 
determined by the above distribution, are given by expressions identical with Eqs. (|8.30[) . 
Then, with both types of initial conditions one obtains the same Langevin equation after 
a time of the order of the width of the memory kernels /C a b (t) , which is the characteristic 
time for the "transient" terms J^b ^ab(t — io)Wb(to) to die out. 
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are sharply peaked at r = 0, the remainder terms in the memory integrals 
change slowly enough in the relevant range, and the kernels enclose a finite 
non-zero algebraic area. Under these conditions, one can replace the kernels 
by Dirac deltas and no memory effects occur. 

Doing this with the memory kernel (18.271) . we write 



/C a b(i~) = 27 ab 5(r) 



(8.31) 



where the 7 a b are damping coefficients related with the strength and charac- 



teristics of the coupling (see below). Then, on using / °°dr 5{r)h{r) 
equation (I8.26P reduces to 



dA 
~dt 



{A,H} + Y,{A,WJ /a(t) + ^7ab 



dWb 
dt 



with 



h(0)/2, 



.32) 



(8.33) 



</a(*)> = , (Ut)f h (t')) = 2 lah k B T5(t - 1') . 

Inspecting Eq. (I8.3ip . one sees that the damping coefficients can be ob- 
tained from the area enclosed by the memory kernels or, alternatively, by 
inserting the definition of the kernel (18.281) into the corresponding integral 
and then using / °°dr cos (cur) = ttS(uj): 



7ab 



dr /C ab (r) 



7ab 



^ tot 



(8.34) 



The area / °°dr /C a b(r) must be: (i) finite and (ii) different from zero, for 
the Markovian approximation to work. The second expression gives the 
damping coefficients in terms of the distribution of normal modes and system- 
environment coupling constants, and could be useful in cases where it could 
be difficult to find the kernels exactly. 



8.4 Examples: Brownian particles and spins 

In order to particularize the general expressions to definite situations, we 
only need to specify the structure of the coupling terms F & . 
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Brownian particle. For instance, let us set F a (p,q) = —c a q (bilinear 
coupling), and write down Eq. (18.241) for A = q and A = p with help from 
{p, B} = —dB/dq and {q, B} = dB/dp. Then one gets dq/dt = dH/dp 
plus Eq. ( 18.31) . which is the celebrated generalized Langevin equation for 
a "Brownian" particle. The fluctuating force is explicitly given by f(t) = 
^2 a c a f a (t) and the memory kernel by /C(t) = ^2 a c 2 JC a {r). Naturally, in 
the Markovian limit /C(t) = 2m / ~f5(r) we have 



dq 
dt 



on 

dp 



dp 
dt 



dU 

dq 



/(t) ~ 7 cE 



^8.35) 



whose relaxation term comprises minus the velocity — {dq/dt) of the particle. 

In general, when {^4, F a } in Eq. ( 18.321) is not constant, the fluctuat- 
ing terms f a (t) enter multiplying the system variables (multiplicative fluc- 
tuations). In this example, owing to the fact that {g, — c a g} = and 
[p, — c a q] = c a , the fluctuations are additive. 



Spin-dynamics. Let us now particularize the above results to the dynam- 
ics of a classical spin. To do so, we merely put A = Si, i = x,y,z, in Eq. 
( I8.24p . and then use Eq. (18. 14j) to calculate the Poisson brackets required. 
Using also dW h /dt = (dW h /ds) ■ (ds/dt), we have 



dsj 
dt 



s A 



&H 

ds 



s A 



ds 



dW h d£ 
~W ' dt 



On gathering these results for % = x, y, z in vectorial form and recalling the 
definition of the effective field B = —&H/ds, we arrive at 



ab 



dt 



ds ds 



Then, defining the fluctuating magnetic field 

dW», 



ds 



(8.36) 



and the second-rank tensor A with elements 

dW a dW h 



A; 



a.b 



7ab- 



ds* dsn 



(8.37) 
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we finally obtain the Langevin equation for the spi 




(8.38) 



Equation (18.381) contains ds/dt on its right-hand side, so it will be referred 
to as a Gilbert-type equation. For e < 1, on replacing perturbatively that 
derivative by its conservative part, ds/dt ~ s AB, one gets the weak-coupling 
Landau-Lifshitz-type equation 



(8.39) 



which describes weakly damped precession. From the statistical properties 
(I8.33P of the fluctuating sources / a (t), one gets 




140) 



which relates the structure of the correlations of the fluctuating field and the 
relaxation tensor. 

For a general form of the spin-environment interaction, due to the oc- 
currence of the tensor A the structure of the relaxation terms in the above 
equations deviates from the forms proposed by Gilbert and Landau and Lif- 
shitz. However, if the spin-environment interaction yields uncorrelated and 
isotropic fluctuations (Ay = X5ij), one finds that: (i) the statistical prop- 



erties (18.401) reduce to those in (18.21) and (ii) the Langevin equation (18.391) 
reduces to the stochastic Landau-Lifshitz equation (18.21) . 

We remark in closing that the occurrence of the vector product s A£ in the 
dynamical equations entails that the fluctuating terms enter in a multiplica- 
tive way. In the spin-dynamics case, in analogy with the results obtained 
for mechanical rigid rotators, the multiplicative character of the fluctuations 
is an inevitable consequence of the Poisson bracket relations for angular- 
momentum-type dynamical variables {sj,^} = J2k e ijk s k, which, even for 
_F a linear in s, lead to non- constant {A, F a } in Eq. flQ4j) . In our derivation 
this can straightly be traced back by virtue of the Hamiltonian formalism 
employed. 



32 Although we omit the symbol of scalar product, the action of a dyadic AB on a 
vector C is the standard one: (AB)C = A(B ■ C). 
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8.5 Discussion 



We have seen that starting from a Hamiltonian description of a classical 
system interacting with the surrounding medium, one can derive generalized 
Langevin equations, which, in the Markovian approach, reduce to known 
phenomenological Langevin equations. 

Note however that the presented derivation of the equations is formal in 
the sense that one must still investigate specific realizations of the system- 
plus-environment whole system, and then prove that the assumptions em- 
ployed (mainly that of Markovian behavior) are at least approximately met. 
Let us give an example for a particle coupled to the elastic waves (phonons) 
of the substrate where it moves. The interaction would be proportional to 
the deformation tensor T~Cse oc du/dx in one dimension. Expanding the 
displacement field in normal modes u{x) = exp(ifcr), where the Uk 

are the coordinates of the environment variables (our Q a ), we have TCse oc 
J2k ikexp(ikx)v,k, so that cx ifc exp(i/cx). If we had allowed complex c a , 
the products c 2 a would had been replaced by \c a \ 2 . Well, the corresponding 
memory kernel [Eq. (18.281) ]. then gives 

/C(r) = e 2 ■ a \ cos(u; q t) Wk —^ f dfc ^ cos(c/cr) oc - D - . 
„ Jo c F r 

But, sin(f2r) / r plays the role of a Dirac delta for any process with time-scales 
much larger than 1/Q. Thus, taking the Markovian limit is well justified in 
this case. 

On the other hand, we have considered the classical regime of the en- 
vironment and the system. A classical description of the environment is 
adequate, for example, for the coupling to low-frequency {hw a /k-oT 1) 
normal modes. Nevertheless, the fully Hamiltonian formalism used, allows 
to guess the structure of the equations in the quantum case (just replacing 
Poisson brackets by commutators). 
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APPENDICES 



Dynamical equations for the averages: macro- 
scopic equation 

From the master equation one can derive the dynamical equations for the av- 
erages of a Markov stochastic process. We shall write down the corresponding 
derivations directly in the multivariate case. 

Let us first write the equation for the time evolution of an arbitrary 
function (/(y))o First, one has 



j t (f(y)) = jJdyf(y)P(y,t) = |dy/(y) 



dy/(y) Jdy'W{y\y')P{y\t) - /dy/(y) Jdy' W(y'\y)P(y,t) 



dP(y,t) 
dt 

Then, by using the master equation to express dPj dt, one has 

|c/(y» = 



,y'<-*y 



r dy7(y) JdyW(y'\y)P(y,t)- jdy f(y) Jdy' W(y'\y)P(y,t) 
= jdyP(y,t)jdy'[f(y')-f(y)}W(y'\y). (.41) 
On applying now this equation to /(y) = yi, and defining [cf. Eq. (14.21) ] 

af\y,t) = Jdy'{ 1 / i - Vi )W(y f \y), (.42) 
the last line in Eq. (j.4ip is the average of a\ (y,t), so we can finally write 



d (y t ) = (a?\y,t)) (i = l,2,...) 



dt 



;.43) 



This is an exact consequence of the master equation and therefore holds for 
any Markov process. 



33 



Here we use the same notation for the stochastic process and its realisations. 
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Note that when is linear function of y one has (a^(y, £)) = a^\(y) , t), 
whence 

which is a system of ordinary differential equations for (y) and can be identi- 
fied with the macroscopic equation of the system. For instance in the decay 
problem (14.91) . since W n ' tTl is non-zero for n' = n — 1, we have 

a^\n,t) = 2j(n - n) W^n = [(^— 1)- i\ jn = -jn . 

Therefore, from the general result (1.431) we have d (n) /dt = (a^{n, t)) = 
— (jn), in agreement with Eq. (14.101) . 

On the other hand, when a\ is a non-linear function of y, Eq. (1.431) is 
not a differential equation for (yi). Then Eq. (I.43p is not a closed equation 
for (y^ but higher-order moments enter as well. Thus, the evolution of (yi) 
in the course of time is not determined by (y^ itself, but is also influenced 
by the fluctuations around this average. To get equations for higher-order 
moments we proceed analogously. For instance, for (yi(t)yj(t)) , we can use 
Eq. (OH) with f(y) = y^. Writting (y-J/J- - !/<%■) = (Ui-ViYiVj ~Vj) +Vi(yj ~ 
Vj) + Vjiv't ~ Vi)i an d defining, analogously to Eqs. (14. 2 p and (l.42p . 

a^(y,t) = J&M-yiWj-VjWtfly). (.44) 

we finally have 



dt 



(.45) 

which is also an exact consequence of the master equation]^ However, if afff 
is a non-linear function of y, the equation involves even higher order moments 



34 Note that for one variable (or for i = j) Eq. (1.45|) reduces to 
±(yi) = (a^(y,t))+2(ya^(y,t)), 

where [cf. Eq. (|3~2)) ] 

a^(y,t)= fdy'(y'-y)W(y'\y) , a^(y,t) = J Ay' (y' - y) 2 W(y'\y) , (.46) 
are the one- variable counterparts of Eqs. (|.42p and (|.44[) . respectively. 
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(yiVjUk), so what we have is an infinite hierarchy of coupled equations for the 
moments. 



The Langevin process £(t) as the derivative of 
the Wiener— Levy process 

Let us formally write dW/dt = £(£), and see which are the properties of the 
W(t) so defined. On integrating over an interval r, we have 

/t + T 
£{s)ds . (.47) 

Let us show that this w(r) is indeed a Wiener-Levy process. Firstly, w{r) is 
Gaussian because £(£) is so. Furthermore, on using the statistical properties 
(JE2D one gets (r,n,r 2 > 0) 

w(0) = , (w(t)) = , (u;(ri)tu(r a )) = 2D minfa, r a ) . (.48) 

Proof: w(0) = follows immediately from the definition (j.47p . while for the 

average (w(r)), one gets (iw(r)) = J t * +T (^(s)) (is = 0. On the other hand, for 

o 

(tu(Ti)tu(T2)), one finds 

2D<5(s-s') 

(w{n)w(T 2 )) = [ t+T1 f +T2 WW^ds'ds 
Jt Jt 

1 

, " » 

/t+min(Ti,T2) /■t+max(ri,r 2 ) 
5{s-s')ds'ds 

= 2Dmin(ri, r 2 ) , 

where we have sorted the integrals to ensure that, when using the Dirac delta 
to take one of them, the location of the "maximum" of the delta is inside the 
corresponding integration interval, and the result is therefore unity. 

Now on comparing these results with those for the increment of the Wie- 
ner-Levy process, whose average is zero since that of the Wiener-Levy pro- 
cess is zero and the second moment is given by Eqs. (I3.12p . one realises that 
the process defined by Eq. ( 1.471) coincides with the increment of a Wiener- 
Levy process^ 

35 They exactly coincide if w(t) is multiplied by \ j\j2D. 
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Proof of the convergence of the Heun scheme 

We shall check that the Heun scheme correctly generates the Kramers-Moyal 
coefficients, by carrying out the Taylor expansion of Eq. (I7.4p . accounting for 
Eq. (17.51) . Concerning the terms involving Ai, one has 

l{Ai(y,t + At)+Ai[y(t),t]}At A,At+j: k B Jk Aw k 

( OA- dA- ' * s 1 

= || A[y(t), t] + -^At + J2d^ " + ' ' ' + My{t),t]jAt 

= A l [y(t),t}At + 0[(Atf 2 }, 
whereas, the terms involving Bi k can be expanded as 
\ {B ik (y, t + At) + B ik [y(t),t)}AW k A d At+j: e B je Aw e 

= i{B ik [y(t),t) + -^At + J2 % ~ VM + ■■■ + B ik [y(t), t]jAW k 

= B ik [y(t),t}AW k + J2^[y(t),t}J2BAy(t),t]AW e AW k + 0[(Atf^ . 

k £ 

In this case we have retained in yj — yj(t) terms up to order (At) 1 / 2 , which 
in the corresponding expansion of Ai are omitted since they yield terms of 
order (At) 3 / 2 . Finally, on inserting these expansions in Eq. (I7.4p . on gets 

Vi {t + At) ~ Vl (t) + Ai\y(t),t]At + B ik[y(t), t]AW k 

k 

+ lj2{T, B dy(t),t}^[y(t),t]\AW k Aw i , (.49) 

kt ^ j ) 

which corresponds to Eq. (2.8) of Ramirez-Piscina, Sancho and Hernandez- 
Machado. Finally, to obtain the Kramers-Moyal coefficients, we have to 
average Eq. (1. 49ft for fixed initial values y(t) (conditional average). To do 
so, one can use (AW*) = and (AW k AW e ) = (2DAt)8 k£ , to get 



J2B ik AW k ) = 0, 
k ' 

jkt yj 1 v jk yj 7 

Y,B lk AW k Y J B j Aw\ = 2D (j2 B ikB jk ) At . 
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Therefore, from Eq. (j.49p one obtains 

( yi (t + At) - w (t)> =^A i + Dj^ B i^§f) At + 0[(Atf 2 ] 

jk ^ 

([ yi (t + At) - Vi {t)] [y 3 (t + At) - Vj (t)}) = 2D(j^B ik B jk ^At + 0[(Atf 2 } 
which lead to the Kramers-Moyal coefficients (15.151) via Eq. (I4.2ip . Q.E.D. 
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Proof of the Box— Muller algorithm. 



We can verify that the transformation (17. 6p leads to a pair of independent 
Gaussian random numbers as an exercise of transformation of variables as 
introduced in Sec. 12.41 



Wi,W 2 



(w 1 ,w 2 ) = dr 1 dr 2 5[wi — \/—2 ln(ri) cos(27rr 2 )] 1 by hypothesis 
Jo Jo , , — 



x 5[w 2 - A/-21n(ri) sm(2vrr 2 )] PR 1 ,R 2 (r 1 ,r 2 ) ■ 
Let us now introduce the substitution 

Ui(r u r 2 ) = A/-21n(ri)cos(27rr 2 ) , u 2 (n,r 2 ) = a/-2 ln(ri) sin(27rr 2 ) , 
the Jacobi matrix of which reads 





dui \ 




1 9n 


dr 2 ' 




I dv.2 


du2 i 




\ dri 


dr2 1 





■± ,\ ( -- cos(27rr 2 ) -2^-2 ln(n) sin(27rr 2 ) 
-^ ^_ 2 1 ln(ri) sin(2 7 rr 2 ) 2^-2 ln(n) cos(27rr 2 ) 

and the corresponding Jacobian (the determinant of this matrix) is given by 
d{u\, u 2 )/d(ri, r 2 ) = —2^jr\. Nevertheless, when changing the variables in 
the above integrals one needs the absolute value of the Jacobian of the inverse 
transformation, which is given by \d(ri,r 2 )/d(ui,u 2 )\ = r\j2-n. Besides, 
ri(ui,u 2 ) can be obtained from the above transformation: — 21n(ri) = u\ + 
u\ =^> r\ = exp[— \ {u\ + w 2 )]. On using all these results the probability 
distribution of (wi,w 2 ) is finally given by 



/oo /»oo i 

dui / du 2 5(w 1 -u 1 )5(w 2 -u 2 ) — exp[-±(uj + ul)] 
■oo </-oo 



' — oo 

1 , , OS 1 



exp (-§w?) -7= exp (-|iWa) • 



2tt v z ' V2tt 

This expression demonstrates that when r x and r 2 are independent random 
numbers uniformly distributed in the interval (0,1), the random variables 
Wi and w 2 given by the transformation (I7.6p are indeed independent and 
Gaussian-distributed with zero mean and variance unity. 



104 



